	********************************************************************************
	**** Date: 5.03.2019														****
	**** Update: 2.09.2022 														****
	**** Author: JW josephgwright@gmail.com 									****	
	** NOTE: This program has been executed in Stata 17.0						****
	***************** **************************************************************

	***********************************************************************************
	**** Additional modules and packages required for executing the following code ****
	***********************************************************************************
		** btscs.ado
		** ciplot.ado
		** estout.ado
		** ivreg2.ado	
		** krls.ado
		** listtex.ado
		** reghdfe.ado
		** sutex.ado
		** xtivreg2.ado
		** crossfold.ado
		** cvauroc.ado

	******************************
	**** Set directory, seed *****
		global dir = "C:\Users\jgw12\Dropbox\Research\Pers-NAVCO\Pers-RENAVCO\CSW-BJPS-reproduction"
		cd "$dir"
		
		capture log close
		log using Pers-Protest.log, replace
		set scheme plotplain
		set more off 
		set matsize 1000
		global seed ="984353"
	
	**********************************************************************
	*********  Gather and merge data sets; create using variables ********
	**********************************************************************
	/*
		* Clean covariate data sets *
		/*
			qui do clean-prio
			qui do clean-nmc
			qui do clean-nelda  /* downloaded 3.8.2018 */
		*/
		* Coup data *
			set more off
			insheet using http://www.uky.edu/~clthyn2/coup_data/powell_thyne_coups_final.txt, clear  /* downloaded 3.8.2018 */
			rename ccode cowcode
			rename country pt_country
			* 1st coup in Yemen (ccode=680) 1968 is South Yemen but the second coup in August is actually in N/All-Yemen
			recode cow (679=678)
			recode cow (680=678) if year==1968 & month==8
			gen  p = "/"
			egen d = concat(day p month p year)
			gen date  = date(d, "DMY")
			gen coupA = coup==1
			gen coupS = coup==2
			drop p d
			local var = "coup coupA coupS"
			foreach v of local var {
				bysort cow year: egen max`v'=max(`v')
				replace `v' =max`v'
				drop max`v'
			}
			egen tag= tag(cow year)
			keep if tag==1
			drop tag
			tsset cow year
			sort cow year
			saveold coups-merge,replace
			
		* EPR data *
			use epr-original, clear
			recode cow (260=255) (679=678) /* Yemen and West Germany are continuous regimes across unification */
			tsset cow year
			gen e = l.ethfrac
			replace ethfrac  = e 
			drop e
			tsset cow year
			gen excluded = l.lrexclpop 
			gen  loggdp = ln(gdpcapl)
			gen logoil = ln(1+oilpcl)
			tsset cow year
			gen gr= d.loggdp
			tssmooth ma  grow = gr, window(2 0 0)
			replace grow =gr if grow==. & gr~=.
			tsset cow year
			gen gr1 = l.gr
			replace gr = gr1
			drop gr1
			sort cow year
			saveold epr-merge,replace
			
		* Counter-weights data * downloaded on 1.25.2019 from:  https://www.ethz.ch/content/dam/ethz/special-interest/gess/cis/international-relations-dam/Data/Coup-Proofing%201970-2017.dta
			use counter-weights-original, clear
			gen cowcode = ccode
			recode cow (679=678) if year>1990
			sort cow year
			saveold	 counter-merge,replace
			use SSFD-cy-jul2020,clear
			gen cowcode=ccode
			recode cow (679=678) if year>1990
			drop cw_* gwf_* milregimetype regime polity2 insample
			local var = "version pcount cbcount ha_cbcount affcount counterbalancing politicization ccount ptotal cbtotal afftotal pg intel police troops militia military border  mean_p mean_cb mean_aff militiachange newmilitia pgchange newpg intelchange newintel policechange newpolice troopschange newtroops milchange newmilitary borderchange newborder"
			foreach v of local var {
				rename `v' debruin_`v'
			}
			sort cow year	
			tab debruin_cbcount
			recode debruin_cbcount (9 8 7 6 5 = 4)
			saveold debruin-merge,replace
		
		* NAVCO data * Version 2.1 (campaign-years from 1945-2013)  /* downloaded 1.27.2020 from https://dataverse.harvard.edu/dataset.xhtml?persistentId=doi:10.7910/DVN/MHOXDV */
			use "NAVCO2-1_ForPublication.dta", clear 
			sum year  /* note that this ends in 2013 */
			gen start = cyear==0 & prim_meth==1
			gen start_regch = cyear==0 & prim_meth==1 & camp_goal==0
			gen ongoing =  cyear>0 & prim_meth==1
			gen violent = 1 if prim_meth==0
			gen nonvio = 1 if prim_meth==1
			gen regch = 1 if camp_goal==0
			gen regchnv = camp_goal==0 & nonvio==1
			gen regchv = camp_goal==0 & violent==1
			gen duration = end_year-start_year+1
			gen cowcode = lccode
			recode cow (347=345) (364=365)        /* Russia/USSR get same cowcode; Yugoslavia */
			recode cow (531=530) if year<1990     /* Eritrea coded as Ethiopia prior to 1990 */
			local var = "duration start start_regch ongoing violent nonvio regch regchnv regchv"
			foreach v of local var {
				egen nav21_`v' = max(`v'), by(cow year)
			}
			egen tag = tag(cow year)
			rename location nav21_country
			keep if tag==1
			keep cowcode year nav21_*
			keep if year<=2010 & year>=1946
			sort cow year
			merge cow year using gwf-original
			tab _merge        /* note that GWF consider pre-1955 Vietnam as "not independent" */
			local var = "duration start start_regch ongoing violent regchnv nonvio regchv"
			foreach v of local var {
				tsset cowcode year
				gen l1_nav21_`v' = l.nav21_`v' 
				recode l1_nav21_`v' (.=0)
				recode nav21_`v' (.=0)
			}		
			drop if _merge==1 
			drop _merge
			tsset cowcode year
			sort cowcode year
			order cowcode year nav21* l1_nav*
			save navco21-merge,replace
			
			
		* NAVCO 1.3 * Version 1.3 (campaigns from 1900-2019)  https://dataverse.harvard.edu/dataset.xhtml?persistentId=doi:10.7910/DVN/ON9XND on 9.21.2020
			import excel using "NAVCO13List.xlsx", clear firstrow		
			gen duration = EYEAR- BYEAR+1
			list CAMPAIGN LOCATION BYEAR EYEAR FAILURE ONGOING   if ONGOING==1  & EYEAR<2018
			drop ONGOING
			expand duration
			sort NAVCOID
			gen n=_n
			egen min=min(n),by(NAVCOID)
			gen year = BYEAR if n==min
			bysort NAVCOID: replace year = year[_n-1] + 1 if n>min & year==.
			gen country  =LOCATION
			replace country="Congo Brazzaville" if country=="Congo-Brazzaville (ROC)"
			replace country="Congo, Dem. Rep." if country=="Democratic Republic of Congo"
			replace country="Germany East" if country=="East Germany"
			replace country="Yemen, Republic of" if country=="Republic of Yemen"
			replace country="Congo Brazzaville" if country=="Republic of the Congo"
			replace country="Yemen Arab Rep" if country=="Yemen Arab Republic"
			replace country="Yemen, (PDR)" if country=="Yemen People's Republic"
			replace country="Congo/Zaire" if country=="Zaire/DRC"
			gen cowcode=.
			replace cowcode =626 if country=="South Sudan"
			qui do cowcodes
			tab LOCATION if cowcode==. & year>=1946
			drop if year<1946 | cowcode==. | year==.
			gen violent = VIOL==1
			gen violent_regch = VIOL==1 & REGCHANGE==1
			gen nonviolent = NONVIOL==1
			gen nonviolent_regch = NONVIOL==1 & REGCHANGE==1
			gen start = nonviolent==1 & year==BYEAR
			gen start_regch = nonviolent==1 & year==BYEAR & REGCHANGE==1
			gen ongoing = nonviolent==1 & year>BYEAR & year<=EYEAR	
			gen p="-"
			egen idyear = concat(cowcode p year)
			local var = "duration start start_regch ongoing violent violent_regch nonviolent nonviolent_regch"
			foreach v of local var {
				egen nav13_`v' = max(`v'),by(idyear)
			}
			egen tag=tag(idyear)
			gen nav13_country = country
			sum cowcode year nav13_*
			keep if tag==1
			keep cowcode year nav13_*
			keep if year<=2010 & year>=1946
			tsset cowcode year
			sort cowcode year
			merge cow year using gwf-original
			tab _merge
			local var = "duration start start_regch ongoing violent violent_regch nonviolent nonviolent_regch"
			foreach v of local var {
				tsset cowcode year
				gen l1_nav13_`v' = l.nav13_`v' 
				recode l1_nav13_`v' (.=0)
				recode nav13_`v' (.=0)
			}			
			drop if _merge==1 
			drop _merge
			tsset cowcode year
			sort cowcode year
			order cowcode year nav13* l1_nav*
			save navco13-merge,replace
			
			
		* MEC data *  https://github.com/ulfelder/nonviolent-uprisings-replication/raw/8646b4410983dd4877c32b7bbdeac91d40d3413d/data/nvc.transformed.csv
			use nvc-original, clear
			keep year country sftgcode yrborn yrdied nvcstart nvc* yth* pol* wditrade wdipopurbmi
			gen cowcode =.
			qui do cowcodes
			replace cowcode = 490 if country=="Congo-Kinshasa"
			replace cowcode = 265 if country=="East Germany"
			replace cowcode = 345 if country=="Federal Republic of Yugoslavia"  /* Serbia 1992-2002 */
			replace cowcode = 816 if country=="North Vietnam"
			sort cow year
			gen repeat=1 if cow==cow[_n-1] & year==year[_n-1]
			drop if repeat==1
			xtset cow year
			gen lag_nvcongoing = l.nvcongoing
			sort cow year
			save mec-merge,replace
			 
		* RENAVCO region data *
			use renavco, clear
			gen cowcode =.
			gen country = statenme
			qui do cowcodes
			sum cowcode  ccode
			replace cowcode = 490 if country=="Democratic Republic of the Congo"
			replace cowcode = 265 if country=="German Democratic Republic"
			replace cowcode = 817 if country=="Republic of Vietnam"
			replace cowcode = 678 if country=="Yemen Arab Republic"
			replace cowcode = 680 if country=="Yemen People's Republic"
			tab country if cowcode==.  /* none are GWF autocracies, 1946-2010 */
			drop if cowcode==.
			sort cow year
			gen repeat = cow==cow[_n-1] & year==year[_n-1]
			list country year if repeat==1
			drop if repeat==1
			drop repeat
			xtset cow year
			gen lnregion = l.renavco_nvc_onset_lregcnt
			replace lnregion = renavco_nvc_onset_lregcnt if year==1946
			tab lnregion
			gen lnregion2  = lnregion
			replace lnregion2= 1.61 if lnregion>1.6 & lnregion~=.
			keep cowcode year country lnregion renavco_nvc_onset_lregcnt
			sort cow year
			save renavco-region-merge,replace
 
		* Fariss data *   /* downloaded 3.8.2018 */
			use fariss-original,clear
			recode cow (679=678) if year>1990
			local var="ciri disap kill polpris tort amnesty state hathaway itt genocide rummel massive_repression executions killing additive"
			foreach v of local var {
				replace `v' = "." if `v'=="NA"
				destring `v', replace
			}
			sort cow year
			xtset cow year
			gen  r1 = latentmean*-1   /* flip scale */
			tssmooth ma lag_repress=r1, window(3 0 0)
			tssmooth ma lag2_repress=r1, window(2 0 0)
			gen l1r = l1.r1
			gen l2r = l2.r1
			gen l3r = l3.r1
			replace l3r=l2r if l3r==.
			replace l3r=l1r if l3r==.
			replace l2r=l1r if l2r==.
			replace lag_repress =lag2_repress if lag_repress==.
			replace lag_repress=l1r if lag_repress==.
			xtset cow year
			gen lag1repress=l1.r1
			gen lag2repress=l2.r1
			gen lag3repress=l3.r1
			drop r1 lag2_repress
			saveold fariss-merge,replace
				
		* VDem data * /* downloaded 6.26.2020 from https://www.v-dem.net/en/data/data-version-10/ */
			cd "$dir"
			use V-Dem-CY-Full+Others-v10,clear
			keep country_name year COWcode v2x_polyarchy v2x_libdem v2x_partipdem v2x_freexp_altinf v2x_frassoc_thick ///
				v2x_clpol v2x_jucon v2juhcind v2cltort v2clkill v2x_corr e_migdppc e_migdppcln v2caconmob v2cademmob v2caautmob
			rename country_name vdem_country
			rename COWcode cowcode
			recode cow (679=678) if year>1990
			tab vdem_country if cow==.
			drop if cow==.
			tab vdem_country if cowcode==99999   /* cases in our sample */
			drop if cowcode==99999
			tsset cow year
			sort cowcode year
			save "vdem-merge.dta", replace
			use V-Dem-CY-Core-v8_merge,clear
			gen cowcode =ccode
			recode cow (679=678) if year>1990
			sort cowcode year
			merge cow year using vdem-merge
			tab _merge
			rename _merge v2merge
			sum year cow
			xtset cow year
			alpha v2clkill v2cltort, item std gen(vkill)
			replace vkill = vkill*-1
			tssmooth ma lag_vkill= vkill,window(3 0 0)
			local var = "v2caconmob v2cademmob v2caautmob v2x_polyarchy v2x_partipdem v2x_freexp_altinf v2x_frassoc_thick v2x_clpol v2x_jucon v2juhcind e_v2x_neopat"
			foreach v of local var {
				gen l1`v'=l1.`v'
				gen l2`v'=l2.`v'
				gen l3`v'=l3.`v'
				gen l4`v'=l4.`v'
			}
			sort cow year
			save vdem-merge, replace
*/
		* Personalism data *
			use gwf-original,clear			
			* Generate binary variables *
			gen milmerit_persA = milmerit_pers
			recode milmerit_persA (2=1) (1=1)
			gen milmerit_persB = milmerit_pers
			recode milmerit_persB (2=1) (1=0)  
			tsset gwf_caseid year
			gen newparty =support==1 & l.support==0
			tsset gwf_caseid year
			replace newparty=1 if l.newparty==1 & l.gwf_leaderid==gwf_leaderid & year==year[_n-1]+1
			gen createparty =militparty_new==1 | (newparty==1  & partyhistory_post==1) 
			gen partyorigin = ldr_group_domparty==1 | ldr_group_priordem==1
			gen offpersmil = officepers*ldr_group_military
			gen offpersparty = officepers*partyorigin 
			tab offpers*
			
			* Label variables *
			global pvars1 = "partyexcom_pers partyrbr offpersparty createparty"
			global pvars2 = "milnotrial milmerit_persB paramil_pers sectyapp_pers offpersmil"
			label var offpersparty "Appointments to high office, party"
			label var offpersmil "Appointments to high office, military"
			label var createparty "Create new party"
			label var partyexcom_pers "Party exec committee"
			label var partyrbr "Rubber stamp party"
			label var milmerit_persB "Military promotions"
			label var milnotrial "Military purge"
			label var sectyapp_pers "Security apparatus"
			label var paramil_pers "Paramilitary"
			
			* Personalism variables *
			set seed 2453456
			irt (2pl partyexcom_pers partyrbr officepers createparty milnotrial milmerit_persB paramil_pers sectyapp_pers) 
			estat report partyexcom_pers partyrbr officepers createparty milnotrial milmerit_persB paramil_pers sectyapp_pers, byparm sort(b)
			predict pers_2pl, latent se(pers_se_2pl)
			qui sum pers_2pl
			gen xpers = (pers_2pl+abs(r(min))) / (r(max) - r(min))
			hist xpers, bin(50)
			irt (2pl milnotrial paramil_pers sectyapp_pers offpersmil) (grm milmerit_pers)
			predict irtpers2, latent  
			qui sum irt
			replace irt = (irt+abs(r(min))) / (r(max) - r(min))
			alpha $pvars2, std gen(ypers2) item
			qui sum ypers2
			replace ypers2 = (ypers2+abs(r(min))) / (r(max) - r(min))
			alpha $pvars1, std gen(xpers1) item
 			qui sum xpers1
			replace xpers1 = (xpers1+abs(r(min))) / (r(max) - r(min))
			alpha partyexcom_pers partyrbr officepers  createparty pleb, std   item gen(ypers1)
			qui sum ypers1
			replace ypers1 = (ypers1+abs(r(min))) / (r(max) - r(min))
		
			* Linear link for paramil_pers sectyapp_pers to allow their errors to be correlated *
			gsem (PER->paramil_pers sectyapp_pers,reg var(PER@1)vce(cluster gwf_leaderid) ///
				cov(e.sectyapp_pers*e.paramil_pers)) ///
				(PER->milnotrial offpersmil,logit) (PER-> milmerit_pers,ologit)
			predict xpers2, latent ebmeans
			qui sum xpers2
			replace xpers2 = (xpers2+abs(r(min))) / (r(max) - r(min))
			estat ic				
			 * GSEM w/out covariances *
			qui gsem (PER->paramil_pers sectyapp_pers,reg var(PER@1)vce(cluster gwf_leaderid)) ///
				(PER->milnotrial offpersmil,logit) (PER-> milmerit_pers,ologit)
			estat ic
			
			* Additional covariances with sectyapp_pers *  More complex covariance structures do not improve fit 
					* First look at the full Gaussian model with 1 covariance specified, for comparison *
			qui gsem (PER-> sectyapp_pers paramil_pers milmerit_pers milnotrial offpersmil,reg var(PER@1)vce(cluster gwf_leaderid) ///
				cov(e.sectyapp_pers*e.paramil_pers)) 
			estat ic
			local var "milmerit_pers milnotrial offpersmil"
			foreach v of local var {
				qui gsem (PER-> sectyapp_pers paramil_pers milmerit_pers milnotrial offpersmil,reg var(PER@1)vce(cluster gwf_leaderid) ///
					cov(e.sectyapp_pers*e.paramil_pers)cov(e.sectyapp_pers*e.`v')) 
				estat ic
			}
				
			sum xpers1 xpers2 xpers	
			spearman xpers1 xpers2 xpers	
			twoway lpolyci xpers1 xpers2,xtit(Security personalism)ytit(Party personalism)legend(off)
			
			twoway (scatter irtpers xpers2,xtit("Generalized SEM")ytit("Hybrid 2PL + GRM")) (lpoly irtpers xpers2,lcol(blue)lpat(solid)) ///
				(lfit   irtpers xpers2,lpat(solid)lcol(red)legend(lab(2 "nonlinear fit") lab(3 "linear fit") order(2 3) ///
				pos(5)ring(0))text(.7 .3 "{&rho}=0.9915"))
			graph export "$dir\sfp-logit-vs-general.pdf", as(pdf)   replace

			 
			* IRT plots *
			qui irt 2pl partyexcom_pers partyrbr officepers createparty milnotrial milmerit_persB paramil_pers sectyapp_pers
			irtgraph iif    (sectyapp_pers,lcolor(blue)) (milmerit_persB,lcolor(red)) (milnotrial,lcolor(green)) ///
			(paramil_pers,lcolor(cyan)),legend(col(2) pos(6)) title(Security & military items) saving(t2,replace) ///
			ylab(,glcolor(gs15)) xtitle("Personalism ({&theta})")
			irtgraph iif  (officepers,lcolor(blue)) (partyexcom_pers,lcolor(red)) (partyrbr,lcolor(green)) ///
			(createparty,lcolor(cyan)),legend(col(2) pos(6))  title(Party & personnel items) saving(t1,replace) ///
			ylab(,glcolor(gs15)) xtitle("Personalism ({&theta})")
			gr combine t1.gph t2.gph, col(2)   ysize(5.5) xsize(9)  ycommon
			erase t1.gph
			erase t2.gph

			* Variance decomposition *
			qui xtset gwf_caseid year
			foreach v in xpers xpers1 xpers2 ypers2 irtpers {
				qui xtsum  `v'  
				qui scalar sdb = r(sd_b)
				qui scalar sdw = r(sd_w)
				qui scalar vart= sdb + sdw
				qui scalar varr = sdw / vart
				scalar list sdw
				scalar list varr
			}
			sort cow year
			save temp-pers,replace
			
		* Merge data sets *
			use temp-pers,clear
			merge cow year using fariss-merge
			tab _merge
			drop if _merge==2
			tab year if _merge==1  /* note that 1946-1948 are not in the Fariss data */
			rename _merge merge1
			sort cow year
			save temp,replace
		    merge cow year using navco13-merge
			tab _merge
			drop if _merge==2
 			rename _merge merge2
			sort cow year
			save temp,replace	
			merge cow year using navco21-merge
			tab _merge
			drop if _merge==2
 			rename _merge merge3
			sort cow year
			save temp,replace	
			merge cow year using epr-merge
			tab _merge
			tab gwf_country if _merge==1  /* note that Singapore and Oman are not in EPR data */
			drop if _merge==2
 			rename _merge merge4
			sort cow year
			save temp,replace
			merge cow year using coups-merge
			tab _merge
			tab year if _merge==2  /* note that coup data begin in 1950 */
			drop if _merge==2
 			rename _merge merge5
			recode coup* (.=0) if year>=1950 & year<=2010
			sort cow year
			save temp,replace
			merge cow year using nelda-merge
			tab _merge
			drop if _merge==2
 			rename _merge merge6
			rename nelda_country Nelda_country
			recode nelda* (.=0) if year>=1945 & year<=2010
			rename Nelda_country nelda_country
			sort cow year
			save temp,replace
			merge cow year using prio-mergeB
			tab _merge
			drop if _merge==2
 			rename _merge merge7
			rename prio_country Prio_country
			recode prio* (.=0) if year>=1945 & year<=2010
			rename Prio_country prio_country
			sort cow year
			save temp,replace
			merge cow year using nmc-merge
			tab _merge
			drop if _merge==2
 			rename _merge merge9
			sort cow year
			save temp,replace
			merge cow year using GWFtscs
			tab _merge
			rename _merge merge10
			sort cow year
			merge cow year using latent-protest-data
			tab _merge
			rename _merge merge11
			sort cow year
			merge cow year using nvc-merge
			tab _merge   
			tab gwf_country if _merge==1  
			rename _merge merge12
			tab year if merge12==1
			tab gwf_country if merge12==1 & year==1990 /* 1990 for South Yemen and East Germany not in MEC data */
			sort cow year   
			merge cow year using counter-merge
			tab _merge  if year>=1970 
			rename _merge merge13
			list gwf_country year if merge13==1 & year>=1970
			keep if gwf_caseid~=.
			sort cow year
			merge cow year using renavco-region-merge
			tab _merge   
			rename _merge merge14
			tab gwf_casename if merge14==1
			tab year if merge12==1 & cowcode==265  /* East Germany 1950--1953 */
			tab year if merge12==1 & gwf_country=="Oman"  /* 1946--1970 */
  			keep if gwf_caseid~=.
			sort cow year
			merge cow year using debruin-merge
			tab _merge   
			rename _merge merge15
			tab gwf_casename if merge15==1
  			keep if gwf_caseid~=.
			sort cow year
			merge cow year using vdem-merge
			tab _merge   
			rename _merge merge16
			tab gwf_casename if merge16==1
			drop if merge16==2
  			keep if gwf_caseid~=.
			xtset cow year
			saveold temp,replace
	
		* Variable transformations *
			use temp,clear
			qui sum latentmean 
			gen repression = (latentmean - r(mean))/r(sd)  /*standardize within sample */
			replace repression = repression*-1   /* flip scale */
			gen lt = ln(gwf_leader_duration)
			gen civwar = prio_lconflict_intra>0  & prio_lconflict_int_intra==2
			gen intwar = prio_lconflict_inter>0  & prio_lconflict_int_inter==2
			xtset gwf_caseid year
 			gen coup12 =  l1.coupA==1 | l2.coupA==1  | l1.coupS==1 | l2.coupS==1 if year>=1950
 			gen election = nelda_mpinc_boycott 
			tab gwf_prior
			gen priordem = gwf_prior=="democracy"
			gen institutions = (support*legcomp)/8  /* 0-1 scale */
			
		* Decade dummies *
			gen d1960 = year>=1960 & year<1970
			gen d1970 = year>=1970 & year<1980
			gen d1980 = year>=1980 & year<1990
			gen d1990 = year>=1990 & year<2000
			gen time = year-1946
			gen time2 = time^2
			gen time3 = time^3
			gen d2000 = year>=2000 & year<=2010
			
		* Save data *
			xtset gwf_caseid year
			keep if gwf_caseid~=.
			sort cow year
			saveold temp, replace
		
		* RENAVCO, NAVCO, and MEC data *
			use "$dir\RENAVCO_nonviolent campaigns_crosswalk_gwfautocracies.dta" ,clear
			tab byear if eyear==.
			recode eyear (.=2015)  /* right censored */
			gen spell =  eyear-byear +1
			expand spell
			sort id
			gen n = _n
			egen min =min(n),by(id)
			gen year=byear if min==n
			sort id
			bysort id: replace year = year[_n-1] +1 if year==.
			gen xonset = byear==year
			gen reduration = 1 if xonset==1
			gen cowcode = ccode
			sort id year
			replace reduration = reduration[_n-1]+1 if reduration==.
			egen max = min(reduration), by(cow year)
			replace reduration = max
			egen maxo = max(xonset),by(cow year)
			replace xonset=maxo
			gen lreduration =ln(reduration)
			gen xongoing = 1
			sort cow year
			sort cow year
			gen repeat = year==year[_n-1] & cow==cow[_n-1]
			tab repeat
			egen r = max(repeat),by(cow year)
			egen mm =mean(byear),by(cow year)
 			listtex cow target year campaign if r==1 & byear==year & mm==year using NVC-multiples.tex, rstyle(tabular) ///
				head("\begin{tabular}{l l c l}"" \textit{Country code} & \textit{Target} & \textit{Year} & \textit{Campaign}  \\\\") ///
				foot("\end{tabular}") replace	
				gen c = campaign if xonset==1 & byear==year
			list cow year target campaign byear if r==1 & c=="",clean noobs
			drop if  r==1 & c==""
			egen tag = tag(cow year)
			tab tag
			keep if tag==1		
			drop tag min max* n r mm repeat c
			sort cow year
			save "$dir\renavco-merge.dta", replace
 

	**********************************************************************************************************
			use "$dir\renavco-merge.dta", clear
			merge cow year using "$dir\temp.dta"
			tab _merge
			drop if _merge==1
			drop _merge
			keep if gwf_caseid~=.
			recode xonset (.=0)
			recode xongoing (.=0)
			xtset cow year
			gen lag_xongoing = l.xongoing
			btscs xonset year cow,gen(xyrs)
			gen lxyrs  = ln(1+xyrs)
			gen xyrs2 = xyrs^2 
			gen xyrs3  = xyrs^3
			
			gen gdem = gwf_fail_subsregime==1
			replace gdem =1 if gwf_casename=="Soviet Union 17-91" & year==1991
			replace gdem =1 if gwf_casename=="Germany East 49-90" & year==1990	

			gen armed = gwf_fail_type  
			recode armed (5=1) (6=1) (7=1) (1=0) (2=0) (3=0) (4=0) (8=0) (9=1) /* coups, rebellions, and foriegn invasions, state collapse */
			recode armed (1=0) if (cow==652 & year==1958) | (cow==345 & year==1990)  /* Syria joins Egypt; Yugoslav republics declare independence */
				
 	**************************************
	********* Look at the sample *********
	**************************************
			**** Plot time trend in combined Onsets ***
			egen yr_protest = mean(xonset),by(year)
			twoway (bar yr_protest year, sort bcol(gs10)xlab(1950(10)2010)lcol(gs15)) (lpoly yr_ year,bw(2) ///
				lcol(blue)lpat(solid)ytitle("{it:Share of dictatorships with }" "{it:Non-violent protest campaigns}", ///
				size(small))legend(lab(1 "% per year") lab(2 "Smoothed trend") pos(6)ring(1)col(2)))
 			graph export "$dir\protest-trend.pdf", as(pdf)   replace

			**** Estimating sample ****
			global cvar ="lxyrs lt lnregion xpers1 xpers2 xpers"
			sum xonset 
			reg xonset $cvar  gdem
			gen sample=e(sample)==1
			tab gwf_casename if sample==0
			keep if sample==1
			label var xonset  "Non-viol. protest campaign onset"
			label var xpers "Personalization"
			label var xpers1  "Party personalization" 
			label var xpers2  `""{bf:Security}     " "{bf:personalization}""'
			label var lpopl "Population (log)"
			label var lt  "Leader tenure (log)" 
			label var lnregion `""Region Non-violent " "episode onsets (log)""'  
			label var gwf_pers  "Personalist regime" 
			label var lxyrs "Time since last onset (log)"
			label var coldwar "Cold war"
 
			* Sum stats *
			*sutex xonset year xpers xpers1 xpers2 xyrs lt lpopl lnregion  if sample==1,minmax labels file($dir/SumstatsStandarized.tex) replace
			sort gwf_country year
			replace navco11 ="No" if navco11==""
			replace mec ="No" if navco11==""
			gen xrenavco="Yes" if renavco==1
			replace xrenavco="No" if renavco==0
			*listtex gwf_country year campaign renavco navco11 mec if xonset==1 & sample==1 using nvcstarts.tex, rstyle(tabular) ///
			*	head("\begin{tabular}{rc l c c c}"" \textit{Country} & \textit{Year} & \textit{Campaign} & \textit{RENAVCO} & \textit{NAVCO} & \textit{MEC} \\\\") foot("\end{tabular}") replace
			sum year
			egen tag =tag(gwf_caseid) if sample==1
			egen max =max(year),by(gwf_caseid)
			egen min=min(year),by(gwf_caseid)
			* listtex gwf_casename min max if tag==1 using regimeslist.tex, rstyle(tabular) ///
			*	head("\begin{tabular}{l c c c}"" \textit{Regime-case} & \textit{Begin year} & \textit{End year}") ///
			*	foot("\end{tabular}") replace
			drop tag max min
				
			* Raw data *
			ttest xpers,by(xonset)
			ttest xpers1,by(xonset)
			ttest xpers2,by(xonset)
	
			* Standardize variables *
			local var ="$cvar"
			foreach v of local var {
				qui sum `v' if sample==1
				qui replace `v' = (`v'-r(mean))/r(sd)
			}
		
			* Get unit and time means *
			egen rmean =mean(xonset) if sample==1,by(gwf_caseid)
			gen max = rmean>0 & rmean<1 if rmean~=.
			tab max
			gen fe = max*-1 +1
			tab fe max
			replace fe = gwf_caseid if max==1
			local var = "xonset ypers2 xyrs xyrs2 xyrs3 coldwar ld gdem $cvar "
			foreach v of local var {
				egen mp_`v'=mean(`v') if sample==1,by(fe)
				egen m_`v'=mean(`v') if sample==1,by(gwf_caseid)
				egen y_`v'=mean(`v') if sample==1,by(year)
				gen mXy_`v' = m_`v'*y_`v'
			}	
			
			
 			sort cow year
			save "$dir\temp-fe.dta", replace
			
			
			* First paragraph stats *
			use temp-fe, clear
			tab xonset xongoing  		/* all onsets are ongoing years */
			egen m_xongoing =mean(xongoing),by(gwf_caseid)
			egen tag =tag(gwf_caseid) if m_xongoing~=.
			gen noprotest = m_xongoing==0
			tab noprotest if tag==1  	/* 58% of regimes never face a mass protest */
			gen noonset  = m_xonset==0
			tab noonset if tag==1		/* 60% of regimes never have an onset between 1946 and 2010 */
			sum xongoing   				/* 7.4% of regime-years have a protest */
			
 ************************************************************************************
 ************************************   Analysis   **********************************
 ************************************************************************************
		use temp-fe,clear
		global cvar = "coldwar lxyrs lt lnregion"
		qui centile xpers2 if sample==1,centile(50)
		di r(c_1)
		gen treat  = xpers2>=-.02472032 if xpers2~=.
		tab treat if sample==1
		egen m_treat = mean(treat) if sample==1,by(gwf_caseid)
		egen min=min(year) if treat==1,by(gwf_caseid)
		gen first_treatment =year==min if year~=.
		sort cowcode year
		keep if year~=. & xpers2~=. & gdem~=.
		save temp-id,replace
		
		* Show within equivalency in OLS *
			xtset gwf_caseid year
				* unit means *
			qui reg xonset lxyrs lt lnregion treat coldwar m_lxyrs m_lt m_lnregion m_coldwar m_treat m_xonset, vce(cluster gwf_caseid)
			lincom treat
				* unit intercepts *
			qui reg xonset lxyrs lt lnregion treat coldwar i.gwf_caseid, vce(cluster gwf_caseid)
			lincom treat		
				* fixed effects *
			qui xtreg xonset lxyrs lt lnregion treat coldwar,fe vce(cluster gwf_caseid)
			lincom treat
				* absorbed fixed effects *
			qui reghdfe xonset lxyrs lt lnregion treat coldwar,a(gwf_caseid) vce(cluster gwf_caseid)
			lincom treat
				* unit means without \bar{Y} *
			qui reg xonset lxyrs lt lnregion treat coldwar m_lxyrs m_lt m_lnregion m_coldwar m_treat, vce(cluster gwf_caseid)
			lincom treat
			* Nonlinear link *
			qui probit xonset lxyrs lt lnregion treat coldwar m_lxyrs m_lt m_lnregion m_coldwar m_treat, vce(cluster gwf_caseid)
			margins,dydx(treat)
			lincom treat
			qui xtprobit xonset lxyrs lt lnregion treat coldwar m_lxyrs m_lt m_lnregion m_coldwar m_treat, vce(cluster gwf_caseid)
			margins, dydx(treat) predict(pu0)
			lincom treat
 			
			 * ECM LRM *
			xtset gwf_caseid year
			ivreg  xonset (d.xonset=l.xonset) d.lxyrs lxyrs d.lt lt d.lnregion lnregion d.treat treat, cluster(gwf_leaderid)
			ivreg  xonset (d.xonset=l.xonset) d.lxyrs lxyrs d.lt lt d.lnregion lnregion d.xpers2 xpers2, cluster(gwf_leaderid)
			
			* LDV-FE * 
			xtset cowcode year
			gen l1mean5 = l.mean5
			alpha l1v2cademmob l1mean5 lag_xongoing  l1_nav13_ongoing l1_nav21_ongoing, std item gen(l1protest)
			gen l2protest = l.l1protest
			reghdfe xonset lxyrs l1protest l2protest lt lnregion treat,a(gwf_caseid year)vce(cluster gwf_leaderid)
			reghdfe xonset lxyrs l1protest l2protest lt lnregion xpers2,a(gwf_caseid year)vce(cluster gwf_leaderid)
			
			* FD *
			xtset gwf_leaderid year
			xi:xtivreg2 xonset lxyrs lt lnregion treat,fd cluster(gwf_leaderid)  
			xi:xtivreg2 xonset lxyrs lt lnregion xpers2,fd cluster(gwf_leaderid)  
			
			* Leader-FE  *
			reghdfe xonset lxyrs lt lnregion treat,a(gwf_leaderid year)cluster(gwf_leaderid)
			reghdfe xonset lxyrs lt lnregion xpers2,a(gwf_leaderid year)cluster(gwf_leaderid)
			

		* Base model *
		use temp-fe, clear
		xtset gwf_caseid year
		probit xonset $cvar xpers2,vce(cluster gwf_caseid)
		margins,dydx(xpers2 lt lnregion)nose
		est store b0
		xtprobit xonset $cvar xpers2,vce(cluster gwf_caseid)
		margins,dydx(xpers2 lt lnregion)nose
		est store b1
 		xtprobit xonset $cvar xpers2 m_coldwar m_lxyrs m_lt m_xpers2 m_lnregion,vce(cluster gwf_caseid)
		margins,dydx(xpers2 lt lnregion)nose
 		est store b2
		xtprobit  xonset lxyrs lt lnregion xpers2 m_lxyrs m_lt m_xpers2 m_lnregion y_lxyrs y_lt y_xpers2 y_lnregion,vce(cluster gwf_caseid)
		margins,dydx(xpers2 lt lnregion)nose
		est store b3
		xtprobit xonset lxyrs lt lnregion xpers2 m_lxyrs m_lt m_xpers2 m_lnregion y_lxyrs y_lt y_xpers2 y_lnregion mXy_lxyrs mXy_lt mXy_xpers2 mXy_lnregion if sample==1,vce(cluster gwf_caseid)
		margins,dydx(xpers2 lt lnregion)nose 
		est store b4
		
		* Check with IRT-2PL pers instead of generalized *
		qui sum irtpers if sample==1
		replace irtpers=(irtpers-r(mean))/r(sd)
		sum irtpers
		qui egen m_irtpers =mean(irtpers2) if sample==1,by(gwf_caseid)
		qui egen y_irtpers =mean(irtpers2) if sample==1,by(year)
		xtprobit xonset $cvar irtpers2 m_coldwar m_lxyrs m_lt m_irtpers m_lnregion,vce(cluster gwf_caseid)
		margins,dydx(irtpers)nose
		xtprobit xonset lxyrs lt lnregion irtpers2 m_lxyrs m_lt m_irtpers m_lnregion y_lxyrs y_lt y_irtpers y_lnregion,vce(cluster gwf_caseid)
		margins,dydx(irtpers)nose

		
		* Various measures of personalism *
		xtset gwf_caseid year
		probit xonset $cvar gwf_pers,vce(cluster gwf_caseid)
		est store p1
		probit xonset $cvar xpers,vce(cluster gwf_caseid)
		margins,dydx(xpers lt)
		est store p2
		xtprobit xonset  $cvar xpers,vce(cluster gwf_caseid)
		margins,dydx(xpers lt)
		est store p3
		xtprobit xonset $cvar xpers1,vce(cluster gwf_caseid)
		margins,dydx(xpers1 lt)
		est store p4	
		xtprobit xonset $cvar xpers2,vce(cluster gwf_caseid)
		margins,dydx(xpers2 lt lnregion)
		est store p5			
		xtprobit xonset $cvar xpers1 xpers2,vce(cluster gwf_caseid)
		est store p6
		xtprobit xonset $cvar gwf_pers,vce(cluster gwf_caseid)
		
		* Random intercept *
		xtset gwf_caseid year
		melogit xonset $cvar xpers2||gwf_caseid:,vce(cluster gwf_caseid)	
		* Random slope *
		melogit xonset $cvar xpers2||gwf_caseid:xpers2,vce(cluster gwf_caseid)
		melogit xonset $cvar xpers2||gwf_caseid:lt lxyrs,vce(cluster gwf_caseid)
		
		* Drop controls
		 xtset gwf_caseid year
		 qui xtprobit xonset lxyrs coldwar lt lnregion xpers2,vce(cluster gwf_caseid)	
		 est store c1
		 qui xtprobit xonset lxyrs xpers2,vce(cluster gwf_caseid)
		 lincom xpers2 
		 est store c2
		 qui xtprobit xonset lxyrs lt xpers2,vce(cluster gwf_caseid)
		 lincom xpers2 
		 est store c3
		 qui xtprobit xonset lxyrs lnregion xpers2,vce(cluster gwf_caseid)
		 lincom xpers2		 
		 est store c4
		 qui xtprobit xonset lxyrs coldwar xpers2,vce(cluster gwf_caseid)
		 lincom xpers2
		 est store c5
		 qui xtprobit xonset lxyrs coldwar lnregion xpers2,vce(cluster gwf_caseid)
		 lincom xpers2
		 est store c6
		 qui xtprobit xonset lxyrs coldwar lt lnregion xpers2,vce(cluster gwf_caseid)
		 lincom xpers2
		 est store c7
		 qui xtprobit xonset lxyrs lt lnregion xpers2,vce(cluster gwf_caseid)
		 lincom xpers2
		 est store c8	
		 qui xtprobit xonset lxyrs coldwar xpers2,vce(cluster gwf_caseid)
		 lincom xpers2
		 est store c9
		* No controls estimates table *
		estout c1 c2 c3 c4 c5 c6 c7 c8 c9 using TableB2.tex, cells(b(star  fmt(%9.3f)) se(par fmt(%9.2f))) ///
			stats(r2 N N_clust) style(tex) replace label starlevels(* 0.05) title(\label{tabB2}) 
		 
		**** Adding covariates, one at a time ****
		use temp-fe,clear
		spearman xpers2 v2x_clpol v2clkill v2cltort v2juhcind v2x_jucon v2x_frassoc_thick v2x_freexp_altinf v2x_partipdem v2x_libdem v2x_polyarchy
		gen n=_n
		gen beta=.
		gen hi=.
		gen lo=.
		gen hi90=.
		gen lo90=.
		gen varname=""
		xtset gwf_caseid year
		global cvar ="coldwar lxyrs lt lnregion"
		global varA= "lag_xongoing loggdp logoil support gwf_mil leadermil seizure_coup coup12 election ythbul4 wdipopurbmi wditrade"
		global varB= "polparcomp polpolcomp polexconst v2x_polyarchy v2x_libdem v2x_partipdem v2x_freexp_altinf v2x_frassoc_thick "
		global varC= "v2x_clpol v2x_jucon v2juhcind e_v2x_neopat priordem legcomp excluded monoethnic  multiethnic civwar grow"
		global varD= "v2cltort v2clkill lag_repress leadermil milethnic_homo nmc_logmilper nmc_logmilex effectivenumber debruin_cbcount debruin_ha_cbcount lpop"
		local var = "$varA $varB $varC $varD"
		local i =1
		foreach v of local var {
			di "`v'"
			qui xtset gwf_caseid
			qui xtprobit xonset $cvar `v' xpers2,vce(cluster gwf_caseid)	 
			qui nlcom _b[xpers2],post
			matrix beta =e(b)  
			local b = beta[1,1]
			qui replace beta=`b' if n==`i'
			matrix var = e(V) 
			local se =var[1,1]
			qui replace hi = `b' + sqrt(`se')*1.96 if n==`i'
			qui replace lo = `b' - sqrt(`se')*1.96 if n==`i'
			qui replace hi90 = `b' + sqrt(`se')*1.65 if n==`i'
			qui replace lo90 = `b' - sqrt(`se')*1.65 if n==`i'
			qui replace varname = "`v'" if n==`i'
			local i = `i' +1
		}
		label define varlab 1 "Ongoing protest campaign" 2 "GDP per capita (log)" 3 "Oil per capita (log)"  ///
			4 "Support party" 5 "Military junta" 6 "Military leader" 7 "Coup seizure" 8 "Recent coup" 9 "Election" ///
			10 "Youth bulge" 11 "Urban pop." 12 "Trade"   13 "Parcomp" 14 "Polcomp" 15 "Exec. constr." ///
			16 "V-Polyarchy" 17 "V-Liberal" 18 "V-Participation" 19 "V-Free express. info" ///
			20 "V-Free express. org" 21 "V-Civil lib." ///
			22 "V-Judical constr." 23 "V-Judicial indep." 24 "V-Neopatrimonialism" 25 "Prior democracy" ///
			26 "Leg. comp." 27 "Excluded pop" ///
			28 "Monoethnic party" 29 "Multiethnic party" 30 "Civil war" 31 "Growth"  32 "V-Torture" ///
			33 "V-Kill" 34 "Repression" 35 "Military leader" ///
			36 "Ethnic homo military" 37 "Military personnel (log)" 38 "Military exp. (log)" ///
			39 "Effective number"  40 "Counter-weights" 41 "Heavily armed cw" 42 "Population",replace
		label values n varlab
		twoway (scatter beta n if n<=42,mcol(blue)yline(-.1275468,lcol(red)lpat(dash_dot))) ///
			(rspike hi lo n if n<=42,lw(vthin)lcol(blue)) ///
			(rspike hi90 lo90 n if n<=42,lcol(blue)lw(medium)ytitle("{&beta}{sub:Security personalization}", ///
			size(large)height(-2)) ///
			xtitle(Added covariate)yline(0,lpat(dash)lcol(gs6))xlab(1(1)42,valuelabel angle(90))legend(off))
		graph export "$dir\baseline-covariates.pdf", as(pdf)   replace

		* Calendar time trends *
		xtset gwf_caseid year
		xtprobit xonset lxyrs lt lnregion xpers2 if year<1989,vce(cluster gwf_caseid)	
		est store time1 
		xtprobit xonset lxyrs lt lnregion xpers2 if year>=1989,vce(cluster gwf_caseid)	
		est store time2
		xtprobit xonset lxyrs time* lt lnregion xpers2,vce(cluster gwf_caseid)	
		est store time3
 		xtprobit xonset lxyrs period* lt lnregion xpers2,vce(cluster gwf_caseid)	
		est store time4
 		xtprobit xonset lxyrs d19* d20  lt lnregion xpers2,vce(cluster gwf_caseid)	
		est store time5		
		
		* Logits *
		xtset gwf_caseid year
		logit xonset $cvar xpers2,vce(cluster gwf_caseid)
		est store logit1
		xtlogit xonset $cvar xpers2,vce(cluster gwf_caseid) 
		est store logit2
		replace lpop = lpop*5
		xtlogit xonset $cvar xpers2 m_coldwar m_lxyrs m_lt m_xpers2  m_lnregion,vce(cluster gwf_caseid)
		est store logit3
		margins,dydx(xpers2)
		firthlogit xonset $cvar xpers2 m_coldwar mp_lxyrs mp_lt mp_xpers2 mp_lnregion mp_xonset
		est store logit4
		margins,dydx(xpers2)
		
		* Cook et al approach to unit intercepts placed in a probit *
		 probit xonset i.fe $cvar xpers2,vce(cluster gwf_caseid)
		 
		* Drop each region *
		use temp-fe,clear
		xtset gwf_caseid year
		local region = "cacar casia ceeurope easia meast nafrica samerica ssafrica weu"
		local i =1
		foreach v of local region {
			qui xtprobit xonset $cvar xpers2 if region~="`v'",re vce(cluster gwf_caseid)
			lincom xpers2
			est store region`i'
			local i =`i' +1
		 }
		 
		* Drop each decade *
		  qui xtprobit xonset $cvar xpers2 if year>=1960,re vce(cluster gwf_caseid)
		  lincom xpers2
		  est store decade1
		  local decade = "1960 1970 1980 1990 2000"
		  local i =2
		  foreach v of local decade {
			qui xtprobit xonset $cvar xpers2 if d`v'==0,re vce(cluster gwf_caseid)
			lincom xpers2
			est store decade`i'
			local i =`i' +1
		  }
	  
		  * Duration time *
		  xtprobit xonset $cvar xpers2,vce(cluster gwf_caseid)
		  xtprobit xonset coldwar xyrs* lt lnregion xpers2,vce(cluster gwf_caseid)
		  gen lnXxpers =lxyrs*xpers2
		  qui xtprobit xonset coldwar lxyrs lt lnregion xpers2 lnXxpers,vce(cluster gwf_caseid)
		  test lnXxpers
		  gen y1Xxpers = xpers2*xyrs
		  gen y2Xxpers = xpers2*xyrs2
		  gen y3Xxpers = xpers2*xyrs3
		  qui xtprobit xonset coldwar xyrs* y1X y2X y3X lt lnregion xpers2,vce(cluster gwf_caseid)
		  test y1 y2 y3
		  
		  * Alternative coding of DV *
		  use temp-fe,clear
		  gen xonsetNAV13 =  nav13_start==1 
		  gen xonsetNAV21 =  nav21_start==1 
		  gen xonsetMEC = nvcstart 
		  * This comes from the Ulfelder and Chenoweth 2015 data:  
		  * "Major Episodes of Contention (MEC) data set to identify the onset of thesenonviolent episodes"
		  gen xonsetREN = xonset*renavco 
		  recode xonsetREN (.=0)
		  tab xonset xonsetNAV13
		  tab xonset xonsetNAV21
		  tab xonset xonsetMEC
		  tab xonsetMEC xonsetNAV13
		  tab xonsetMEC xonsetNAV21
		  
		  xtset gwf_caseid year
		  btscs xonsetMEC year cow,gen(yrsMEC)
		  btscs xonsetNAV13 year cow,gen(yrsNAV13)
		  btscs xonsetNAV21 year cow,gen(yrsNAV21)
		  btscs xonsetREN year cow,gen(yrsREN)

		  sum yrs*
		  gen lyrsNAV13 =ln(yrsNAV13+1)
		  gen lyrsNAV21 =ln(yrsNAV21+1)
		  gen lyrsMEC =ln(yrsMEC+1)
		  gen lyrsREN =ln(yrsREN+1)

		  gen yrsMEC2 = yrsMEC^2
		  gen yrsMEC3 = yrsMEC^3
		  gen yrsNAV212 = yrsNAV21^2
		  gen yrsNAV213 = yrsNAV21^3
		  gen yrsNAV132 = yrsNAV13^2
		  gen yrsNAV133 = yrsNAV13^3
		  gen yrsREN2 = yrsREN^2
		  gen yrsREN3 = yrsREN^3

		  qui:xtprobit xonsetMEC lyrsMEC period* lnregion lt lpop xpers2,vce(cluster gwf_caseid)
		  lincom xpers2
		  est store alt1
		  qui:xtprobit xonsetMEC yrsMEC* period* lnregion lt xpers2,vce(cluster gwf_caseid)
		  lincom xpers2
		  qui:xtprobit xonsetNAV13 lyrsNAV13 period* lnregion lt xpers2,vce(cluster gwf_caseid)
		  lincom xpers2
		  est store alt2
		  qui:xtprobit xonsetNAV13 yrsNAV13* period* lnregion lt xpers2,vce(cluster gwf_caseid)
		  lincom xpers2
		  qui:xtprobit xonsetNAV21 lyrsNAV21 period* lnregion lt xpers2,vce(cluster gwf_caseid)
		  lincom xpers2
		  est store alt3
		  qui:xtprobit xonsetNAV21 yrsNAV21* period* lnregion lt xpers2,vce(cluster gwf_caseid)
		  lincom xpers2
		  qui:xtprobit xonsetREN lyrsREN period* lnregion lt xpers2,vce(cluster gwf_caseid)
		  lincom xpers2
		  est store alt4
		  qui:xtprobit xonsetREN yrsREN* period* lnregion lt xpers2,vce(cluster gwf_caseid)
		  lincom xpers2
		  qui:xtprobit xonset  lxyrs period* lnregion lt xpers2,vce(cluster gwf_caseid)
		  est store alt0
		  
		  
		  **********************************
		  ************ Plots ***************
		  **********************************  
		  	* Reported RE & CRE models coefficient plot *
			coefplot (b0, msym(S)) (b1, msym(d)) (b2, msym(t)) (b3, msym(O)) (b4, msym(oh)), ///
				title("Personalization and Non-violent protest", size(small)) ///
				drop(_cons xyr* lxyr* m_* mXy_* y_* ) ///
				order(gwf_personal xpers xpers1 xpers2) xline(0) grid(glcolor(gs13)) ///
				mfcolor(white) xlabel(-.8(.2).4,labsize(vsmall))levels(95 90)  ///
				legend(lab(3 "Probit")lab(6 "RE probit")lab(9 "CRE" "probit") ///
				lab(12 "Two-way CRE" "probit")lab(15 "Interactive CRE"  "probit")  ///
				size(vsmall)pos(7)ring(0)) ysize(3) xsize(3.5) saving(r1, replace) ///
				xtitle("Coefficient estimate", size(small) height(6)) ///
				note("90 (thick) and 95 (thin) percent confidence intervals", size(vsmall) pos(6))
				graph export "$dir\base-models.pdf", as(pdf) replace
		  	
			* Personalism variables models coefficient plot *
			coefplot (p1, msym(d)) (p2, msym(t)) (p3, msym(oh)) (p4, msym(plus)) (p5, msym(P)) (p6, msym(T)), ///
				title("Personalization and Non-violent protest", size(small)) drop(_cons xyr* lxyr* ) ///
				order(gwf_personal xpers xpers1 xpers2) xline(0) grid(glcolor(gs13)) mfcolor(white) xlabel(-.4(.1).3,labsize(vsmall))  levels(95 90)  ///
				legend(lab(3 "Regime type" "logit") lab(6 "Personalization" "logit") lab(9 "Personalization"  "RE logit") lab(12 "Party pers." "RE logit") ///
				lab(15  "Security pers." "RE logit") lab(18 "Both types" "personalism") size(vsmall)) ysize(3) xsize(3.5) saving(r1, replace)	xtitle("Coefficient estimate", size(small) height(6)) ///
				note("90 (thick) and 95 (thin) percent confidence intervals", size(vsmall) pos(6))
				graph export "$dir\personalism-models.pdf", as(pdf)   replace
			
				
			* Logit models coefficient plot *
			coefplot (logit1, msym(d)) (logit2, msym(t)) (logit4, msym(plus)), ///
				title("Personalization and Non-violent protest", size(small)) drop(_cons lxyr* m* coldwar ) ///
				order(xpers2 lt lpopl lnregion) xline(0) grid(glcolor(gs13)) mfcolor(white) xlabel(-.75(.25).75,labsize(vsmall))  levels(95 90)  ///
				legend(lab(3 "Logit") lab(6 "RE" "logit")   lab(9 "Firth" "logit") size(vsmall))ysize(3) xsize(3.5) ///
				saving(r1, replace)	xtitle("Coefficient estimate", size(small) height(6)) note("90 (thick) and 95 (thin) percent confidence intervals", size(vsmall) pos(6))
				graph export "$dir\logit-models.pdf", as(pdf)   replace

				
				* LPMs(WFE from Imai and Kim 2019 uses R script pasted below) *
					 use temp-fe,clear
					 keep if sample==1
					 * set the binary treatment variable at the median of security personalism *
					 qui sum xpers2 if sample==1,detail
					 local m=r(p50)
					 gen treat1=xpers2>=(`m'-.0001) 
					    * Note: using the IRT-2PL instead of mixed yields same treatment variable 
						qui sum irtpers2 if sample==1,detail
						local m=r(p50)
						di `m'
						gen treat2=irtpers2>=(`m'-.0001) 
						tab treat2 treat1
						drop treat2 
					 xtset gwf_caseid year
					 gen l1treat1=l.treat1
					 gen l2treat1=l2.treat1
					 gen l3treat1=l3.treat1

					 gen l1xonset=l.xonset
					 gen l2xonset=l2.xonset
					 gen l3xonset=l3.xonset
					 
					 * Bivariate *
					 ttest xonset,by(treat1)
					 ttest xongoing,by(treat1)

					 * Comparison margins estimate from RE probit *
					 global cvarlpm="lnregion lt lxyrs"
					 xtset gwf_caseid year
					 qui xtprobit xonset period* $cvarlpm treat1,  vce(cluster gwf_caseid)
					 margins,dydx(treat1)
							
					 * 2-way FE *
					 qui reghdfe xonset $cvarlpm treat1,a(gwf_caseid year)vce(cluster gwf_leaderid)
					 lincom treat1
					 est store lpm1

					 * Unit-specific quadratic time trend *
					 set matsize 10000
					 xtset gwf_caseid year
					 xi:qui reg xonset i.gwf_caseid*time i.gwf_caseid*time2 ///
						lxyrs lt lnregion treat1,vce(cluster gwf_leaderid)
					 lincom treat1
					 est store lpm2
					 
					 * Interactive FE *
					 qui regife xonset $cvarlpm treat1,a(gwf_caseid year)factor(gwf_caseid year,1) vce(cluster gwf_leaderid)
					 lincom treat1
					 est store lpm3
					 
					 * Past treatment *
					 reghdfe xonset $cvarlpm treat1 l1treat1 l2treat1,a(gwf_caseid year)cluster(gwf_leaderid)
					 est store lpm4
					 lincom treat1
					 test l1treat1 l2treat1  

					 * Past outcome *
					 reghdfe xonset $cvarlpm treat1 l1xonset l2xonset,a(gwf_caseid year)cluster(gwf_leaderid)
					 est store lpm5
					 lincom treat1
					 test l1xonset l2xonset  
					 
					 * Past outcome, low level *
					 reghdfe xonset $cvarlpm treat1 l1v2cademmob l2v2cademmob,a(gwf_caseid year)cluster(gwf_leaderid)
					 est store lpm6
					 * Check with 4 lags, per Hamilton 2018 RES
					 reghdfe xonset $cvarlpm treat1 l1v2cademmob l2v2cademmob l3v2cademmob,a(gwf_caseid year)cluster(gwf_leaderid)
					 reghdfe xonset $cvarlpm treat1 l1v2cademmob l2v2cademmob l3v2cademmob l4v2cademmob,a(gwf_caseid year)cluster(gwf_leaderid)

					 * Check with continous *
					 reghdfe xonset $cvarlpm xpers2 l1v2cademmob l2v2cademmob l3v2cademmob,a(gwf_caseid year)cluster(gwf_leaderid)
					 reghdfe v2cademmob $cvarlpm treat1 l1v2cademmob l2v2cademmob l3v2cademmob,a(gwf_caseid year)cluster(gwf_leaderid)
					
					 * Past outcome and past treatment *
						* First test whether t-2 and t-3 lags of outcome/treatment are related to outcome conditional on 1-year lags *
					 reghdfe xonset $cvarlpm treat1 l2xonset l3xonset l2treat1 l3treat1,a(gwf_caseid year) cluster(gwf_leaderid)
 					 test l2xonset l3xonset l2treat1 l3treat1
					 reghdfe xonset $cvarlpm treat1 l1treat l1xonset l2xonset l3xonset l2treat1 l3treat1,a(gwf_caseid year) cluster(gwf_leaderid)
					 test l2xonset l3xonset l2treat1 l3treat1
					 ivreghdfe xonset $cvarlpm treat1 (l1treat1 l1xonset= l2xonset l3xonset l2treat1 l3treat1),a(gwf_caseid year) cluster(gwf_leaderid)
					 lincom treat1
					 est store lpm7
					 
					 * Check with continous treatment *
					 ivreghdfe xonset $cvarlpm xpers2 (l1.xpers2 l1xonset= l2xonset l3xonset l2.xpers2 l3.xpers2), ///
						a(gwf_caseid year) cluster(gwf_leaderid)
					 
					* LPM table *
					estout lpm* using TableB1.tex, cells(b(star  fmt(%9.3f)) se(par fmt(%9.2f))) ///
						stats(r2 N N_clust widstat) style(tex) replace label starlevels(* 0.05) title(\label{tabB1})		
					saveold temp-fe1,replace version(12)

						/*  R script for weight-FE (WFE) estimator: 
						ls()
						setwd("C:/Users/jgw12/Dropbox/Research/Pers-NAVCO/Pers-RENAVCO/CSW-BJPS-reproduction")
						getwd()

						rm(list=ls())
						set.seed(23789)

						## install missing packages, and update if newer version available

						packageList<-c("car","plm","xtable","wfe","sandwich","stargazer")

						for(i in 1:length(packageList)){
							if (!require(packageList[i],character.only = TRUE)){
								install.packages(packageList[i], repos="http://lib.stat.cmu.edu/R/CRAN/")
							}
						}
						# update.packages(ask = FALSE, dependencies = c('Suggests'), oldPkgs=packageList, repos="http://lib.stat.cmu.edu/R/CRAN/")
						 
						library(car)
						library(foreign)
						library(xtable)
						library(wfe)
						library(sandwich)
						library(stats4)
						 

						rm(list=ls())
						set.seed(89270258)
						data.temp<-read.dta("temp-fe1.dta")

						# Unweighted unit FE + period effects #
						mod.W1<-wfe(xonset~treat1+lxyrs+lnregion+lt+
									  period1+period2+period3+period4+period5+
									  period6+period7+period8+period9+period10+
									  period11+period12,
									 data=data.temp,treat="treat1",
									 unit.index="gwf_caseid",time.index="year",
									 method="unit", qoi="ate",hetero.se=TRUE,
									 auto.se=TRUE,unweighted=TRUE)
						 summary(mod.W1)
						 
						 # Weighted unit FE + period effects #
						 mod.W2<-wfe(xonset~treat1+lxyrs+lnregion+lt+
									   period1+period2+period3+period4+period5+
									   period6+period7+period8+period9+period10+
									   period11+period12,
									  data=data.temp,treat="treat1",
									  unit.index="gwf_caseid",time.index="year", 
									 method="unit",qoi="ate",hetero.se=TRUE,
									  auto.se=TRUE,unweighted=FALSE)
						summary(mod.W2)
						 

						# Two-way FE UNweighted #
						mod.W3<-wfe(xonset~treat1+lxyrs+lnregion+lt,
									data=data.temp,treat="treat1",
									unit.index="gwf_caseid",
									time.index="year", method="unit",
									qoi="ate",hetero.se=TRUE,
									auto.se=TRUE,unweighted = TRUE,
									estimator = "did")
						summary(mod.W3)

						# Two-way FE Weighted #
						mod.W4<-wfe(xonset~treat1+lxyrs+lnregion+lt,
									data=data.temp,treat="treat1",
									unit.index="gwf_caseid",
									time.index="year", method="unit",
									qoi="ate",hetero.se=TRUE,
									auto.se=TRUE,unweighted =FALSE,
									estimator = "did")
						summary(mod.W4)

						*/
						
			* Semiparametric estimator *
			use temp-fe,clear
			xtset gwf_caseid year
			set seed 897243574
			xi:qui xtsemipar xonset i.year xpers2 lt lxyrs lnregion,nonpar(xpers2)cluster(gwf_leaderid)deg(3)gen(a b)spline
			qui gen invb= invlogit(b) /* rescale on 0,1 */
			qui sum invb if sample==1
			di r(mean)
			qui sum xonset   if sample==1 & invb~=.
			replace invb = invb-(.49943119-r(mean)) if sample==1 
			sum invb xonset if sample==1 & invb~=.
			qui gen inva=invlogit(a)
			qui sum inva if sample==1
			di r(mean)
			qui sum xonset if sample==1 & inva~=.
			replace inva = inva-(.4983158-r(mean)) if sample==1 
			twoway (hist xpers2,bin(50)bcol(gs12)yscale(range(0 3000)axis(2))yaxis(2) ///	
				ylab(none,axis(2))freq ytitle("",axis(2)))(lpolyci invb xpers2,bw(.5) xtit(Security personalism)ylab(.03(.01).06) yscale(alt) ///
				legend(off) ytitle(Probability of campaign onset)col(blue)xlab(-1.5(.5)1.5)yline(.0390294)) 
 			graph export "$dir\onset-semipar.pdf", as(pdf) replace
			drop a b invb inva
			
			* Time trend models coefficient plot *
			coefplot (time1, msym(d)) (time2, msym(t)) (time3, msym(oh)) (time4, msym(plus)) (time5, msym(P)), ///
				title("Adjusting for time trends", size(small)) drop(_cons xyr* lxyr* time* d* coldwar period*) ///
				order(xpers2) xline(0) grid(glcolor(gs13)) mfcolor(white) xlabel(-.5(.25).5,labsize(vsmall))  levels(95 90)  ///
				legend(lab(3 "1946-1988" "period") lab(6 "1989-2010" "period") lab(9 "Five-year"  "periods") lab(12 "Non-linear"  "time trend")  ///
				lab(15 "Decade"  "dummies") size(vsmall)) ysize(3) xsize(3.5) saving(r1, replace)	///
				xtitle("Coefficient estimate", size(small) height(6))note("90 (thick) and 95 (thin) percent confidence intervals", size(vsmall) pos(6))
				graph export "$dir\time-models.pdf", as(pdf) replace
						
			* Drop each region *
			  local region = "cacar casia ceeurope easia meast nafrica samerica ssafrica weu"
			  local i =1
			  foreach v of local region {
				qui xtprobit xonset lxyrs lt lnregion xpers2 if region~="`v'",re vce(cluster gwf_caseid)
				lincom xpers2
				est store region`i'
				local i =`i' +1
			  }
			  
			* Drop each decade *
			  qui xtprobit xonset lxyrs lt lnregion xpers2 if year>=1960,re vce(cluster gwf_caseid)
			  lincom xpers2
			  est store decade1
			  local decade = "1960 1970 1980 1990 2000"
			  local i =2
			  foreach v of local decade {
				qui xtprobit xonset lxyrs lt lnregion xpers2 if d`v'==0,re vce(cluster gwf_caseid)
				lincom xpers2
				est store decade`i'
				local i =`i' +1
			  }
			  
			label var xpers "Personalization"
			label var xpers1 `""Party        " "personalization""'
			label var xpers2 `""Security    " "personalization""'
			label var lpopl "Population (log)"
			label var lt `""Leader" "tenure (log)""'
			label var lnregion `""Region NVC" "onsets (log)""'
			label var gwf_pers `""Personalist" "regime    ""'
			  
			* Drop regions models coefficient plot *
			coefplot (region1, msym(d)) (region2, msym(t)) (region3, msym(oh)) (region4, msym(plus)) (region5, msym(P)) ///
				(region6, msym(d)) (region7, msym(t)) (region8, msym(oh)) (region9, msym(plus)) , ///
				title("Drop one region at a time", size(small)) drop(_cons pyr* time* d* coldwar) ///
				order(xpers2) xline(0) grid(glcolor(gs13)) mfcolor(white) xlabel(-.3(.15).3,labsize(vsmall))  levels(95 90)  ///
				legend(lab(3 "Central" "America") lab(6 "Central" "Asia") lab(9 "Central/East"  "Europe") lab(12 "East"  "Asia") ///
				lab(15 "Middle"  "East") lab(18 "North"  "Africa") lab(21 "South"  "America") lab(24 "Sub-Saharan"  "Africa") ///
				lab(27 "Western"  "Europe") size(vsmall)) ysize(3) xsize(3.5) saving(r1, replace)	///
				xtitle("Coefficient estimate", size(small) height(6))note("90 (thick) and 95 (thin) percent confidence intervals", ///
				size(vsmall) pos(6))
				graph export "$dir\drop-region-models.pdf", as(pdf) replace
				
			* Drop decade models coefficient plot *
			coefplot (decade1, msym(d)) (decade2, msym(t)) (decade3, msym(oh)) (decade4, msym(plus)) (decade5, msym(P)) (decade6, msym(d)), ///
				title("Drop each decade at a time", size(small)) drop(_cons pyr*) ///
				order(xpers2) xline(0) grid(glcolor(gs13)) mfcolor(white) xlabel(-.3(.15).3,labsize(vsmall))  levels(95 90)  ///
				legend(lab(3 "1940s" "1950s") lab(6 "1960s") lab(9 "1970s") lab(12 "1980s") ///
				lab(15 "1990s") lab(18 "2000s") size(vsmall)) ysize(3) xsize(3.5) saving(r1, replace) xtitle("Coefficient estimate", ///
				size(small) height(6))note("90 (thick) and 95 (thin) percent confidence intervals", size(vsmall) pos(6))
				graph export "$dir\drop-decade-models.pdf", as(pdf) replace
				
			* Alternative DV models coefficient plot *
			coefplot (alt0, msym(d)) (alt1, msym(t)) (alt2, msym(oh)) (alt3, msym(plus)) (alt4, msym(S)), ///
				title("Security personalization and Non-violent protest" "campaign onsets, alternative DV tests", ///
				size(small)) drop(_cons xyrs* yrs* lyrs* period*) order(xpers2 lt lpopl lnregion) ///
				xline(0) grid(glcolor(gs13)) mfcolor(white) xlabel(-.3(.15).3,labsize(vsmall)) levels(95 90) ///
				legend(lab(3 "Original") lab(6 "MEC") lab(9 "NAVCO 1.3") lab(12 "NAVCO 2.1") lab(15 "NEVER") size(vsmall)) ///
				ysize(3) xsize(3.5) xtitle("Coefficient estimate", size(small) height(6))  ///
				note("90 (thick) and 95 (thin) percent confidence intervals", size(vsmall) pos(6))
				graph export "$dir\alt-models.pdf", as(pdf)   replace
				
		*** Cox model with shared frailty ***
			use temp-fe,clear
			gen oxyrs = xyrs+1
			sum oxyrs
			stset oxyrs, fail(xonset)
			stcox coldwar lt lnregion xpers2 lpop,shared(gwf_caseid) nohr
			estat phtest,detail
			gen lnoxyrs = ln(oxyrs)
			gen stXlt = lnoxyrs*lt
			stcox lt coldwar lnregion xpers2 lpop stXlt,shared(gwf_caseid) nohr
			estat phtest,detail	
			centile xpers2,centile(10 90)
			centile lnregion,centile(10 90)
			stcurve, hazard at1(xpers2=-1.6) at2(xpers2=1.3) title(Loyalty) ylab(0(.1).4)  ///
				xtitle(Time since last protest onset,height(6)) saving(r1.gph,replace) ///
				legend(lab(1 "Low loyalty") lab(2 "High loyalty") pos(5) col(1) ring(0)) range(1 40)  
			stcurve, hazard  at1(lnregion=-.7) at2(lnregion=1.8) title(Regional protest) ylab(0(.1).4) ///
				xtitle(Time since last protest onset,height(6)) saving(r2.gph,replace) ///
				legend(lab(1 "Low regional protest") lab(2 "High regional protest") pos(5) col(1) ring(0)) range(1 40) 
			gr combine r1.gph r2.gph,  
			graph export "$dir\cox-models.pdf", as(pdf)   replace
			
		*** IVSLS test ***
			use temp-fe,clear
			** Leader effects **
			xtset gwf_caseid year
			qui xtprobit xonset coldwar lxyrs lt lnregion xpers2, re vce(cluster gwf_caseid)
			margins,dydx(xpers2)
			xtset gwf_leaderid year
			qui xtlogit xonset coldwar lxyrs lt lnregion xpers2, re vce(cluster gwf_leaderid)
			margins,dydx(xpers2)
			qui xtlogit xonset coldwar lxyrs lt lnregion xpers2, fe  
			margins,dydx(xpers2)
				* Leader-specific time trend *
			xi:qui reg xonset i.gwf_leaderid*time lxyrs lt lnregion xpers2,cluster(gwf_leaderid)
			lincom xpers
			
			** IV **
			egen rcount = count(year),by(gwf_caseid)
			keep if rcount>=2
			xtset gwf_caseid year
			gen milit2 = militrank^2
			xtreg xonset i.year lxyrs lt lnregion xpers2,re vce(cluster gwf_caseid)
			est store iv1
			xtivreg xonset i.year lxyrs lt lnregion (xpers2=militrank milit2), ///
				re vce(cluster gwf_caseid)  ec2sls 
			est store iv2
				qui: xtreg xpers2 i.year lxyrs lt lnregion   militrank milit2,vce(cluster gwf_caseid)
				test militrank=0
				test milit2=0,ac
			reghdfe xonset lxyrs lt lnregion xpers2,a(gwf_caseid year) cluster(gwf_leaderid) 
			est store iv3
			xi:ivreg2h xonset i.year lxyrs lt lnregion (xpers2=militrank milit2), ///
				fe cluster(gwf_leaderid) gmm2s partial(i.year)
			est store iv5
			
			 * Within IV-probit *
			ivprobit xonset $cvar (xpers2=militrank milit2) m_coldwar m_lxyrs  m_lt m_xpers2 m_lnregion,vce(cluster gwf_leaderid)
			margins,dydx(xpers2)predict(pr)
			
			 ** IV tests table **
 			estout iv1 iv2 iv3 iv5 using TableB3.tex, cells(b(star  fmt(%9.4f)) se(par fmt(%9.4f))) ///
				stats(N N_clust widstat) style(tex) replace label starlevels(* 0.05) title(\label{tabB3})
		
		*******************************
		*** RMSE model fit analysis ***
		*******************************
			use temp-fe,clear
			xtset cowcode year
			gen l1poldurable = l.poldurable
			gen l1v2x_libdem=l.v2x_libdem
			gen l1ythbul = l.ythbul
			gen l1v2cltort=l1.v2cltort
			gen l1v2clkill=l1.v2clkill
			recode debruin_affcount (5 4 3=3)
			global seed = 2892742
			global k=10 /* folds Chenoweth and Ufelder use 5 folds*/
			global y="xonset"
			global varA= "lag_xongoing l1v2cademmob loggdp logoil support lt leadermil coup12 election ythbul4 wdipopurbmi wditrade"
			global varB= "polparcomp polpolcomp polexconst l1v2x_polyarchy l1v2x_libdem l1v2x_partipdem l1v2x_freexp_altinf l1v2x_frassoc_thick "
			global varC= "l1v2x_clpol l1v2x_jucon l1v2juhcind e_v2x_neopat lnregion legcomp excluded monoethnic  multiethnic civwar grow"
			global varD= "l1v2cltort l1v2clkill lag_repress l1poldurable milethnic_homo lag_milpersonnel lag_milspend effectivenumber debruin_cbcount debruin_ha_cbcount lpop"
			local var = "$varA $varB $varC $varD"
			sum `var'
			qui reg $y xyrs*
			keep if e(sample)==1  
			sutex year $y  $varA $varB $varC $varD xpers2, minmax    	
				xtset gwf_caseid year
			capture program drop crossvalrmse
			program define crossvalrmse
				mat A=r(est)
				mat U = J(rowsof(A),1,1)
				mat c = U'*A
				mat cm = c/rowsof(A)
				mat list cm
			end
			
			gen base=.
			gen rmse=.
			gen n=_n
			gen varname = ""
			local var = " $varA $varB $varC $varD xpers2"	
			local i =1
			foreach v of local var {
				di "`v'"
				qui set seed $seed
				qui crossfold reghdfe $y xyrs* if `v'~=.,a(gwf_caseid coldwar)  k($k) 
				crossvalrmse
				local r = cm[1,1]
				qui replace base = `r'   if n==`i'
				qui set seed $seed
				qui crossfold reghdfe $y xyrs* `v',a(gwf_caseid coldwar)  k($k) 
				crossvalrmse
				local r = cm[1,1]
				qui replace rmse = `r'   if n==`i'
				qui replace varname = "`v'" if n==`i'
				local i = `i' +1
			}
		 
			* Change in RMSE from baseline *
			gen ch = (((rmse)-(base))/(rmse))*100
			* Plot change in RMSE *
			sort ch
			gen xn=_n
			replace n = xn
			drop xn
			list n varname ch if n<=43,clean noobs
		
				
				*** Condition of baseline specification covariates ***
				global varA= "lag_xongoing l1v2cademmob loggdp logoil support   leadermil coup12 election ythbul4 wdipopurbmi wditrade"
				global varB= "polparcomp polpolcomp polexconst l1v2x_polyarchy l1v2x_libdem l1v2x_partipdem l1v2x_freexp_altinf l1v2x_frassoc_thick "
				global varC= "l1v2x_clpol l1v2x_jucon l1v2juhcind e_v2x_neopat   legcomp excluded monoethnic  multiethnic civwar grow"
				global varD= "l1v2cltort l1v2clkill lag_repress l1poldurable milethnic_homo lag_milpersonnel lag_milspend effectivenumber debruin_cbcount debruin_ha_cbcount lpop"	
				local var = "$varA $varB $varC $varD xpers2"	
				local i =1
				foreach v of local var {
					qui set seed $seed
					qui crossfold reghdfe $y xyrs* lnregion lt if `v'~=.,a(gwf_caseid coldwar)  k($k) 
					qui crossvalrmse
					local r = cm[1,1]
					qui replace base = `r'   if n==`i'
					qui set seed $seed
					qui crossfold reghdfe $y xyrs* lnregion lt `v',a(gwf_caseid coldwar)  k($k) 
					qui crossvalrmse
					local r = cm[1,1]
					qui replace rmse = `r'   if n==`i'
					qui replace varname = "`v'" if n==`i'
					local i = `i' +1
				}
				replace ch = (((rmse)-(base))/(rmse))*100
				sort ch
				replace n =_n
				label define vlab 1 "Election" 2  "{bf:Security personalism}" 3 "Monoethnic party"  ///
					4 "Legislative competitiveness" 5 "Multiethnic party" 6 "V-Civil liberties" 7 "V-Free expression info" ///
					8 "Civil conflict" 9 "V-Free expression orgs" 10 "V-Polyarchy"	11 "Log GDPpc" 12 "Support party" ///
					13 "V-Participation"  14 "Ethnically homogenous military" 15 "V-Judical constraint"  ///
					16 "V-Kill" 17 "V-Judicial independence" 18 "Ongoing protest campaign" 19 "V-Liberal democracy" ///
					20 "V-Torture" 21 "Oil per capita (log)" 22 "Military leader" 23 "Excluded pop" 24 "V-Neopatrimonialism" ///
					25 "Recent coup" 26 "Economic growth" 27 "State repression" 28 "Military spending (log)" ///
					29 "Military personnel (log)" 30"Population" 31 "Parcomp" 32 "Polcomp"  33 "Executive constraint"   ///
					34 "Prior pro-democracy mobilization" 35 "Polity durable" 36 "Urban population" ///
					37 "Youth buldge"  38 "Trade"  39 "Counter-weights" 40 "Heavily armed counter-weights" ///
					41 "Effective # military organizations",replace
				label values n vlab,
				list varname ch if n<=41,clean noobs
				twoway (scatter n ch if n<=41,ysize(10)msym(P)mcol(blue)xtitle("Change in RMSE") ///
					ytitle("",height(0)) xlab(0 (2) 10)  ///
					ylab(1(1)41,valuelabel)legend(off)xline(0,lcol(red)))
				graph export "$dir\rsme-onset.pdf", as(pdf)   replace
				
			*********************************************************
			*** Cross-Validation with 10 runs of 5-folds 		  ***
			*** for 41 variables added separately to the baseline ***
			*********************************************************	
				use temp-fe,clear
				drop groups
					gen l1poldurable = l.poldurable
					gen l1polparcomp = l.polparcomp
					gen l1polpolcomp = l.polpolcomp
					gen l1polexconst = l.polexconst
					gen l1v2x_libdem=l.v2x_libdem
					gen l1ythbul4 = l.ythbul4
					gen l1v2cltort=l1.v2cltort
					gen l1v2clkill=l1.v2clkill
					gen l1wdipopurbmi = l1.wdipopurbmi
					gen l1wditrade = l1.wditrade
					gen l1loggdp = l1.loggdp
					gen l1lpop = l1.lpop
				global s = 97454
				global y="xonset"
				global x="lxyrs lt lnregion coldwar"
				global k="5"
				global varA= "xpers2 l1v2cademmob l1loggdp l1lpop election l1ythbul4 lag_milpersonnel lag_milspend l1wdipopurbmi l1wditrade "
				global varB= "lag_xongoing l1polparcomp l1polpolcomp l1polexconst l1v2x_polyarchy l1v2x_libdem"
				global varC= "l1v2x_partipdem l1v2x_freexp_altinf l1v2x_frassoc_thick l1v2x_clpol l1v2x_jucon l1v2juhcind"
				global varD= "e_v2x_neopat legcomp excluded monoethnic multiethnic civwar grow l1v2cltort l1v2clkill lag_repress"
				global varE= "l1poldurable milethnic_homo  effectivenumber debruin_cbcount debruin_ha_cbcount logoil support leadermil coup12"
				
				gen auc0=.
				gen auc1=.
				gen ratio=.
				gen varname=""
				gen n=_n
				local var = "$varA $varB $varC $varD $varE"	
				local i =1
				foreach v of local var {
				    di "`v'"
					qui cvauroc $y $x if `v'~=.,k($k)seed($s)
					local r = r(mean_auc)
					qui replace auc0 = `r' if n==`i'
					qui cvauroc $y $x `v' if `v'~=.,k($k)seed($s)
					local r = r(mean_auc)
					qui replace auc1 = `r' if n==`i'
					qui replace varname = "`v'" if n==`i'
					qui xtset gwf_caseid year
					qui xtsum `v'
					qui replace ratio = r(sd_w)/(r(sd_w)+r(sd_b)) if n==`i'
					local i = `i' +1
				}
				gen ch = (((auc1)-(auc0))/(auc0))*100
				sort ch
				replace n =_n
				browse n varname auc0 auc1 ch 
				twoway (scatter ch ratio) (lpoly ch ratio)
				
				replace varname ="{bf:Security Personalism}" if varname=="xpers2"
				replace varname ="State repression" if varname=="lag_repress"
				replace varname ="Election" if varname=="nelda_election"
				replace varname ="Polity durable" if varname=="l1poldurable"
				replace varname ="V-Polyarchy" if varname=="l1v2x_polyarchy"
				replace varname ="V-Civil liberties" if varname=="l1v2x_clpol"
				replace varname ="Intnl conflict" if varname=="prio_lconflict_int_inter"
				replace varname ="Civil conflict" if varname=="prio_lconflict_int_intra"
				replace varname ="Ethnic military" if varname=="milethnic_homo"
				replace varname ="GDPpc" if varname=="l1loggdp"
				replace varname ="Counter balance" if varname=="debruin_affcount"
				replace varname ="# Security orgs" if varname=="effectivenumber"
				replace varname ="Youth bulge" if varname=="l1ythbul4"
				replace varname ="Population" if varname=="l1lpop"
				replace varname ="Urban population" if varname=="l1wdipopurbmi"
				replace varname ="Support party" if varname=="support"
				replace varname ="Excluded pop" if varname=="excluded"
				replace varname ="Prior coup" if varname=="coup12"
				replace varname ="Oil per capita (log)" if varname=="logoil"
				replace varname ="Legislative competitiveness" if varname=="legcomp"
				replace varname = "Multiethnic party" if varname=="multiethnic"
				replace varname = "Monoethnic party" if varname=="monoethnic"
				replace varname = "Trade" if varname=="l1wditrade"
				replace varname = "Military personnel" if varname=="lag_milpersonnel"
				replace varname = "Military spending" if varname=="lag_milspend"
				replace varname = "Prior pro-democracy mobilization" if varname=="l1v2cademmob"
				replace varname = "Election" if varname=="election"
				replace varname = "V-Free expression info" if varname=="l1v2x_freexp_altinf"
				replace varname = "Military leader" if varname=="leadermil"
				replace varname = "V-Free expression orgs" if varname=="l1v2x_frassoc_thick"
				replace varname = "Counter-weights" if varname=="debruin_cbcount"
				replace varname = "Heavily armed counter-weights" if varname=="debruin_ha_cbcount"
				replace varname = "Economic growth" if varname=="grow"
				replace varname = "Polity parcomp" if varname=="l1polparcomp"
				replace varname = "Polity polcomp" if varname=="l1polpolcomp"
				replace varname = "Polity exconst" if varname=="l1polexconst"
				replace varname = "V-Neopatrimonialism" if varname=="e_v2x_neopat"
				replace varname = "V-Participation" if varname=="l1v2x_partipdem"
				replace varname = "Prior ongoing protest" if varname=="lag_xongoing"
				replace varname = "Civil war" if varname=="civwar"
				replace varname = "V-Judicial independence" if varname=="l1v2juhcind"
				replace varname = "V-Judical constraint" if varname=="l1v2x_jucon"
				replace varname = "State torture" if varname=="l1v2cltort"
				replace varname = "State killing" if varname=="l1v2clkill"
				replace varname = "V-Liberal democracy" if varname=="l1v2x_libdem"		 
				twoway (scatter n ch if n>40 & n<=41,mlab(varname)mlabpos(left)msym(P)mcol(blue)) ///
					(scatter n ch if n<=40,mlab(varname)xscale(range(-4 12.5)) ///
					msym(P)mcol(blue)xlab(-4(4)12)legend(off)ylab(1(1)41,  ///
					labsize(tiny)valuelabel)ytitle(Added variables,size(small)) xline(0,lcol(red))  ///
					xtitle(Percent change in baseline test AUC,size(small)) title("Protest campaign onset")   ///
					text(6 9 "{it:Baseline AUC: 0.614}", size(small)) ysize(8))
				graph export "$dir\crossval-onset.pdf", as(pdf)   replace
	
				* Show that military variables don't increase AUC substantially when population is in the specification *
				local var ="lag_milsp lag_milper"
				foreach v of local var {
					qui cvauroc $y $x lpop if `v'~=.,k($k)seed($s)
					local r1 = r(mean_auc)
					qui cvauroc $y $x `v' lpop if `v'~=.,k($k)seed($s)
					local r2 = r(mean_auc)
					local ch = (`r2'-`r1')/(`r1')
					local s ="  "
					di "`v'" 
					di `r1'  
					di `r2' 
					di `ch'
				}
			
			

				
	******************************************************************
	**************************  REPRESSION  **************************
	******************************************************************
		***** Differenced repression model *****
 		use temp-fe,clear
		xtset gwf_caseid year
		gen ongXpers=xongoing*xpers2
		gen onsXpers = xonset*xpers2
		xi:reghdfe d.repress d.xongoing d.xpers2,absorb(period* lt) cluster(gwf_caseid)
		est store rd1
		xi:reghdfe d.repress d.xongoing d.xpers2 d.ongXpers,absorb(period* lt) cluster(gwf_caseid)
		est store rd2
		centile xpers2 if e(sample)==1,centile(10 90)
		lincom d.xongoing + d.ongX*-1.6
		lincom d.xongoing + d.ongX*1.3
		xi:reghdfe d.repress d.xongoing d.coupA d.coupS d.election d.civwar d.xpers2 lxyrs,absorb(period* lt) cluster(gwf_caseid)
		est store rd3
		xi:reghdfe d.repress d.xongoing d.coupA d.coupS d.election d.civwar d.xpers2 d.ongXpers lxyrs,absorb(period* lt) cluster(gwf_caseid)
		est store rd4
		lincom d.xongoing + d.ongX*-1.6
		lincom d.xongoing + d.ongX*1.3

			** Check split sample **
			xi:qui reghdfe d.repress d.xongoing d.coupA d.coupS d.election d.civwar lxyrs if xpers2~=.,absorb(period* lt) cluster(gwf_caseid)
			lincom d.xongoing
			qui sum xpers2 if e(sample)==1,detail
			local m = r(p50)
			xi:qui reghdfe d.repress d.xongoing d.coupA d.coupS d.election d.civwar lxyrs if xpers2>=`m',absorb(period* lt) cluster(gwf_caseid)
			lincom d.xongoing	
			est store rd5
			xi:qui reghdfe d.repress d.xongoing d.coupA d.coupS d.election d.civwar lxyrs if xpers2<`m',absorb(period* lt) cluster(gwf_caseid)
			lincom d.xongoing
			est store rd6
			
		* Differenced estimates table *
		estout rd1 rd2 rd3 rd4 rd5 rd6 using TableC1.tex, cells(b(star  fmt(%9.4f)) se(par fmt(%9.3f))) ///
			stats(r2 N N_clust) style(tex) replace label starlevels(* 0.05) title(\label{tabC1})
			
		* Additional FEs *
		xi:reg d.repress d.xongoing d.coupA d.coupS d.election d.civwar d.xpers2 d.ongXpers lt, cluster(gwf_caseid)
		xi:reghdfe d.repress d.xongoing d.coupA d.coupS d.election d.civwar d.xpers2 d.ongXpers lt,absorb(year cow) cluster(gwf_caseid)
		xi:reghdfe d.repress d.xongoing d.coupA d.coupS d.election d.civwar d.xpers2 d.ongXpers lt,absorb(year gwf_caseid) cluster(gwf_caseid)
		xtset gwf_leaderid year
		xi:reghdfe d.repress d.xongoing d.coupA d.coupS d.election d.civwar d.xpers2 d.ongXpers lt,absorb(year) cluster(gwf_caseid)
		lincom d.xongoing + d.ongX*-1.6
		lincom d.xongoing + d.ongX*1.3

 		***** Repression in onset year & during protest campaigns *****
 		use temp-fe,clear
		xtset gwf_caseid
		egen mean_repress=mean(repress),by(gwf_caseid)
		egen mx_repress=mean(repress) if xongoing==0,by(gwf_caseid)
		egen m_repress=max(mx_),by(gwf_caseid)
		global x1 = "lpop lnregion lt lxyrs"
		centile xpers2 if xonset==1,centile(50)
		centile xpers2,centile(50)

		qui sum xpers2 if lpopl~=.,detail
		local m=r(p50)
		gen treat1=xpers2>=`m' 		 
		 ciplot lag1repress if xonset==1,by(treat) xtitle(Security personalization,size(small)height(3)) ///
			 xlab(2  "Low"  5  "High"  ) note("")  ///
			 ytitle("Lagged repression",height(2)) title("Repression{sub:t-1}") saving(r1.gph,replace)
		 ciplot lag2repress if xonset==1,by(treat) xtitle(Security personalization,size(small)height(3)) ///
			 xlab(2  "Low"  5  "High"  ) note("")  ///
			 ytitle("Lagged repression",height(2)) title("Repression{sub:t-2}") saving(r2.gph,replace)
		 ciplot lag3repress if xonset==1,by(treat) xtitle(Security personalization,size(small)height(3)) ///
			 xlab(2  "Low"  5  "High"  ) note("")  ///
			 ytitle("Lagged repression",height(2)) title("Repression{sub:t-3}") saving(r3.gph,replace)
		 ciplot lag_repress if xonset==1,by(treat) xtitle(Security personalization,size(small)height(3)) ///
			 xlab(2  "Low"  5  "High"  ) note("")  ///
			 ytitle("Lagged repression",height(2)) title("Repression{sub:t-1:t-3}") saving(r4.gph,replace)
		gr combine r1.gph r2.gph r3.gph r4.gph, ycommon col(2) iscale(.7) ysize(5) xsize(5)
		graph export "$dir\repression-ttests.pdf", as(pdf)   replace

		gen beta = .
		gen n =_n
		gen hi = .
		gen lo=.
		
		krls repress xpers2 $x1 period* lag1repress lag2repress lag3repress if xonset==1  
			mat o=e(Output)
			replace beta = o[1,1] if n==3
			replace hi = o[1,1]+1.95*o[1,2]  if n==3
			replace lo = o[1,1]-1.95*o[1,2]  if n==3
		krls repress xpers2 $x1 period* lag1repress lag2repress if xonset==1
			mat o=e(Output)
			replace beta = o[1,1] if n==2
			replace hi = o[1,1]+1.95*o[1,2]  if n==2
			replace lo = o[1,1]-1.95*o[1,2]  if n==2
		krls repress xpers2 $x1 period* lag1repress if xonset==1
			mat o=e(Output)
			replace beta = o[1,1] if n==1
			replace hi = o[1,1]+1.95*o[1,2]  if n==1
			replace lo = o[1,1]-1.95*o[1,2]  if n==1
		krls repress xpers2 $x1 period* mean_repress if xonset==1 /*condition on average regime repression level in all years */
			mat o=e(Output)
			replace beta = o[1,1] if n==4
			replace hi = o[1,1]+1.95*o[1,2]  if n==4
			replace lo = o[1,1]-1.95*o[1,2]  if n==4
		krls repress xpers2 $x1 period* m_repress if xonset==1    /*condition on average regime repression level in non-NVC years */
			mat o=e(Output)
			replace beta = o[1,1] if n==5
			replace hi = o[1,1]+1.95*o[1,2]  if n==5
			replace lo = o[1,1]-1.95*o[1,2]  if n==5
		krls repress xpers2 $x1 period* if xonset==1    			/*no condition on lagged repression */
			mat o=e(Output)
			replace beta = o[1,1] if n==6
			replace hi = o[1,1]+1.95*o[1,2]  if n==6
			replace lo = o[1,1]-1.95*o[1,2]  if n==6
		browse beta hi lo n if n<=6
		label define varlab 1 "t-1" 2 "t-1,t-2" 3 "t-1,t-2,t-3" 4 "Regime mean" 5 "All non-protest yrs" 6 "No lag",replace
		label values n varlab
		twoway (rspike hi lo n if n<=6,lcol(gs1)) (scatter beta n if n<=6,mcol(gs1)msym(O)yline(0,lcol(red)) yline(.0393769,lcol(blue)) ///
			ytitle("{&beta}{sub:Security personalization}",size(large))xtitle(Repression lag,height(3)) ///
			xscale(range(.75 6.25))yscale(range(0 .1))ylab(0(.05).2) xlab(1(1)6,valuelabel angle(45))legend(off))	
		graph export "$dir\repression-lags.pdf", as(pdf) replace
	
		**** Reported Models ****
		xtset cowcode year
		gen lag4repress=l.lag3repress
		krls repress period* lag1repress lag2repress lag3repress xpers2  if xonset==1, 
				est store rep1
		krls repress period* lag1repress lag2repress lag3repress xpers2 $x1 if xonset==1, 
				est store rep2
		krls repress period* lag1repress lag2repress lag3repress xpers2 lreduration $x1  if xongoing==1, d(dx) 
				est store rep3
		krls repress period* lag1repress lag2repress lag3repress lag4repress  xpers2 $x1 if xonset==1, 
		estout rep1 rep2 rep3 using Table1.tex, cells(b(star  fmt(%9.4f)) se(par fmt(%9.3f))) ///
			stats(r2 N N_clust) style(tex) replace label starlevels(* 0.05) title(\label{tab1})
		twoway (hist lreduration,col(gs14)yscale(range(0 20)axis(2))yaxis(2)bin(30)ylab(0(0)0,axis(2)) ytitle("",axis(2))) ///
			(lpolyci dx_xpers2 lreduration,bw(.75) lcol(blue*1.2)lpat(solid)col(blue*.25)legend(off) ///
			xtitle("Protest campaign duration (years, log scale)", ///
			size(small)) ytitle(Marginal effect of security personalism) yline(0,lcol(red))yscale(alt)yscale(alt axis(2)) ///
			xlab(0 "1" .69 "2"  1.099 "3"   1.61 "5" 2.079 "8" 2.48 "12")) 	 
		drop dx_* beta n hi lo
		graph export "$dir\pers-repression.pdf", as(pdf) replace
		
		* Check with OLS RE *
		xtset gwf_caseid year
		reg repress lag1repress lag2repress lag3repress period* xpers2 $x1 if xonset==1,cluster(gwf_caseid)  /* 6.3 percent marginal */
		xtreg repress lag1repress lag2repress lag3repress period* xpers2 $x1 if xonset==1,cluster(gwf_caseid)  /* 5.5 percent marginal */
		krls repress lag1repress lag2repress lag3repress period* xpers2 $x1 if xonset==1,   /* 4.8 percent marginal */
		
		* Additional covariates *
		gen beta = .
		gen n =_n
		gen var =""
		gen hi = .
		gen lo=.
		global vn = 25
		recode debruin_ha_cbcount (9 8 7 6 5 4 3 =3)
		local c="lnmembers territorial lag_xongoing civwar election coup coupA coupS coup12 excluded ethfrac milethnic_homo milethnic_het logoil loggdp xpers1 support nmc_logmilex nmc_logmilper  debruin_counterbalancing debruin_cbcount debruin_ha_cbcount effective v2juhcind v2x_jucon"
		local i = 1
		foreach v of local c {
			di "`v'"
			qui xi:krls repress xpers2 period* lag1repress lag2repress lag3repress  $x1 `v' if xonset==1
			mat o=e(Output)
			replace beta = o[1,1] if n==`i'
			replace var = "`v'" if n==`i'
			replace hi = o[1,1]+1.95*o[1,2]  if n==`i'
			replace lo = o[1,1]-1.95*o[1,2]  if n==`i'
			local i  = `i'+1
		}
		label define varlab 1 "Protest size" 2 `""Territorial" "campaign""' 3 `""Ongoing NV" "campaign""' 4 "Civil war" 5 "Election" ///
			6 "Coup attempt" 7 "Failed coup" 8 "Successful coup" 9 "Lag coup counts" 10 "Excluded pop." 11 "Ethnic fraction."   ///
			12 `""Ethnically" "homo military""' 13 `""Ethnically" "hetero military""' 14 "Oil per capita (log)"  15 "GDP per capita (log)"  ///
			16 "Party personalism" 17 "Support party" 18 "Mil. spending" 19 "Mil. personnel" 20  "Counterbalancing"  ///
			21  `""Counterbalance" "count""' 22  `""Counterbalance" "Heavily armed""' 23 `""Effective #" "mil. orgs""' ///
			24  "Judicial indep." 25  "Judicial constraint",replace
		label values n varlab
		twoway (rspike hi lo n if n<=$vn) (scatter beta n if n<=$vn,msym(O)mcol(gs1)yline(0,lcol(red)) yline(.0287,lcol(blue)) ///
			ytitle("{&beta}{sub:Security personalization}",size(large))xtitle(Added variable,height(-6)) ///
			xscale(range(.75 $vn.25))yscale(range(0 .08))ylab(0(.02).08) xlab(1(1)$vn,valuelabel angle(90))legend(off)xsize(8)ysize(4))		
		graph export "$dir\added-variables-repression.pdf", as(pdf) replace
		qui: reg repress lag_repress period* xpers2 $x1 if xonset==1,cluster(gwf_caseid)
		lincom xpers2
		qui: xtreg repress lag_repress period* xpers2 $x1 if xonset==1,cluster(gwf_caseid)
		lincom xpers2

		  * Check with VDem data on repression *
		  spearman vkill repress 
		  spearman vkill repress if xonset==1
		krls vkill period* lag_vkill xpers2 if xonset==1 
				est store rep1a
		krls vkill period* lag_vkill xpers2 $x1 if xonset==1
				est store rep2a
		krls vkill period* lag_vkill xpers2 lreduration $x1 if xongoing==1,d(dx)
				est store rep3a
		twoway lpolyci dx_xpers2 lreduration,bw(.7) legend(off) xtitle("Protest campaign duration (log)", ///
			size(small)) ytitle(Marginal effect of security personalism) yline(0,lcol(red))
		drop dx_*
		estout rep1a rep2a rep3a using TableC2.tex, cells(b(star  fmt(%9.4f)) se(par fmt(%9.3f))) ///
		stats(r2 N N_clust) style(tex) replace label starlevels(* 0.05) title(\label{tabC2})
		
		* Check with additive *
		xtset cow year
		gen ladd= l.additive
		krls additive period* ladd xpers2 $x1 if xonset==1

		**** Test Repression model with other forms of personalization *****  adjust Table 1, column 2 model 
		krls repress period* lag1repress lag2repress lag3repress gwf_pers $x1 if xonset==1
		est store p1
		krls repress period* lag1repress lag2repress lag3repress xpers $x1 if xonset==1
		est store p2
		krls repress period* lag1repress lag2repress lag3repress xpers1 $x1 if xonset==1
		est store p3
		krls repress period* lag1repress lag2repress lag3repress xpers2 $x1 if xonset==1
		est store p4
		krls repress period* lag1repress lag2repress lag3repress xpers1 xpers2 $x1 if xonset==1
		est store p5
		estout p1 p2 p3 p4 p5  using TableC3.tex, cells(b(star  fmt(%9.4f)) se(par fmt(%9.3f))) ///
		stats(r2 N N_clust) style(tex) replace label starlevels(* 0.05) title(\label{tabC3})
					
		*** Repression ECMs ***
		use temp-fe,clear
		gen ongXpers=xongoing*xpers2
		gen onsXpers = xonset*xpers2
		gen nvtime =lreduration
		recode nvtime (.=0)
		local var  = "repress xongoing ongXpers xpers2 civwar election coup coupS coupA lnregion lpop"
		foreach v of local var {
			xtset gwf_caseid year
			qui gen d`v'=d.`v'
			qui gen l`v'=l.`v'
		}
		
		* ECMs *
		global x= "dcivwar lcivwar delection lelection dcoupS lcoupS dcoupA lcoupA"
		reghdfe drepress lrepress dxongoing lxongoing dxpers2 lxpers2 lt lxyrs, ///
			absorb(period* nvtime) vce(cluster gwf_caseid nvtime)
		centile xpers2 if e(sample)==1,centile(50)
		est store r1
		reghdfe drepress lrepress dxongoing lxongoing dxpers2 lxpers2 lt lxyrs $x, ///
			absorb(period* nvtime) vce(cluster gwf_caseid nvtime)
		est store r2
		reghdfe drepress lrepress dxongoing dongXpers lxongoing longXpers dxpers2 lxpers2 lt lxyrs, ///
			absorb(period* nvtime) vce(cluster gwf_caseid nvtime)
		est store r3
		reghdfe drepress lrepress dxongoing dongXpers lxongoing longXpers dxpers2 lxpers2 lt lxyrs $x, ///
			absorb(period* nvtime) vce(cluster gwf_caseid nvtime)
		est store r4
		
				* Check ECM split sample *
				qui reghdfe drepress lrepress dxongoing lxongoing  lt lxyrs if xpers2>.0260124,  absorb(period* nvtime) vce(cluster gwf_caseid)
				lincom dxongoing
				qui reghdfe drepress lrepress dxongoing lxongoing  lt lxyrs if xpers2< .0260124,  absorb(period* nvtime) vce(cluster gwf_caseid)
				lincom dxongoing
				qui reghdfe drepress lrepress dxongoing lxongoing  lt lxyrs $x if xpers2> .0260124,  absorb(period* nvtime) vce(cluster gwf_caseid)
				lincom dxongoing
				qui reghdfe drepress lrepress dxongoing lxongoing  lt lxyrs $x if xpers2<.0260124,  absorb(period* nvtime) vce(cluster gwf_caseid)
				lincom dxongoing
				
		 * ECM plots *
			interflex drepress dxongoing xpers2 lrepress lxongoing lnregion lpop lt lxyrs, ///
				fe(period* nvtime) cluster(gwf_caseid) type(kernel)  bw(5)
			mat e=r(margeff)
			gen g=_n
			gen b=.
			gen se=.
			gen npers=.
 			forval i = 1(1)50 {
  				qui replace npers=e[`i',1] if g==`i'
  				qui replace b=e[`i',2] if g==`i'
  				qui replace se=e[`i',3] if g==`i'
			}
			gen hi5=.
			gen hi10=.
			gen lo5=.
			gen lo10=.
			replace hi5  = b+1.96*se
			replace hi10 = b+1.65*se
			replace lo5  = b-1.96*se
			replace lo10 = b-1.65*se
			graph twoway (hist xpers2,bin(50)freq bcolor(gs14)yscale(off)) ///
				(rspike hi5 lo5 npers if g<=50,lcolor(blue*.35)lwidth(thin)yaxis(2)yscale(alt axis(2)) ///
				xtitle("Security personalization",size(small)height(4)) ///
				yline(0,lp(dash)axis(2)) ytitle("Marginal effect of Protest campaign on Repression", ///
				size(small) height(4)axis(2))ylab(-.12(.04).12,axis(2)) ///
				title("Short-run",size(med))) ///
				(rspike hi10 lo10 npers if g<=50,lcolor(blue*.35)lwidth(medthick)yaxis(2)yscale(alt axis(2))) ///
				(line b npers if g<=50,yaxis(2)yscale(alt axis(2))lpattern(solid)lcolor(blue*1.1) ///
				xscale(range(0 1)) xlab(-1.6(.8)1.6)legend(lab(1 "Security personalization") ///
				lab(2 "95% CI") lab(4 "Marginal effect") size(small)pos(6)col(3)ring(1)) ///
				legend(order(1 2 4))saving(g1.gph,replace)) 
			drop se b g npers 
			interflex drepress lxongoing xpers2 lrepress dxongoing lnregion lpop lt lxyrs, ///
				fe(period* nvtime) cluster(gwf_caseid) type(kernel)  bw(5)
			mat e=r(margeff)
			gen g=_n
			gen b=.
			gen se=.
			gen npers=.
 			forval i = 1(1)50 {
  				qui replace npers=e[`i',1] if g==`i'
  				qui replace b=e[`i',2] if g==`i'
  				qui replace se=e[`i',3] if g==`i'
			}
			replace hi5  = b+1.96*se
			replace hi10 = b+1.65*se
			replace lo5  = b-1.96*se
			replace lo10 = b-1.65*se
			graph twoway (hist xpers2,bin(50)freq bcolor(gs14)yscale(off)) ///
				(rspike hi5 lo5 npers if g<=50,lcolor(blue*.35)lwidth(thin)yaxis(2)yscale(alt axis(2)) ///
				xtitle("Security personalization",size(small)height(4)) ///
				yline(0,lp(dash)axis(2)) ytitle("Marginal effect of Protest campaign on Repression", ///
				size(small) height(4)axis(2))ylab(-.12(.04).12,axis(2)) ///
				title("Long-run",size(med))) ///
				(rspike hi10 lo10 npers if g<=50,lcolor(blue*.35)lwidth(medthick)yaxis(2)yscale(alt axis(2))) ///
				(line b npers if g<=50,yaxis(2)yscale(alt axis(2))lpattern(solid)lcolor(blue*1.1) ///
				xscale(range(0 1)) xlab(-1.6(.8)1.6)legend(lab(1 "Security personalization") ///
				lab(2 "95% CI") lab(4 "Marginal effect") size(small)pos(6)col(3)ring(1)) ///
				legend(order(1 2 4))saving(g2.gph,replace)) 
		gr combine g1.gph g2.gph,xsize(8) 
		graph export "$dir\interflex-repression.pdf", as(pdf) replace  
		drop se b g npers 

		*** Onset year relative to post-onset years ***
		xi:krls repress xpers2 xonset civwar lpop lnmembers lt if xongoing==1,d(dx) 
		twoway lpolyci dx_xonset xpers2,legend(off) yline(0,lcol(red))
		drop dx*
			
		******* NAVCO repression **********
		use temp-fe,clear
		gen hipers = xpers2>0 if xpers2~=.
		tab any_repression if xonset==1

		 * Onset years *
		tab any_repression hipers if xonset==1,col
		krls any_repression xpers2 if xonset==1  
		krls any_repression xpers2 lnregion lpop lnmembers lt period* if xonset==1,d(k1)  
		krls any_repression xpers2 lnregion lpop lnmembers lt period* lag1repress lag2repress lag3repress if xonset==1
		twoway lpolyci k1_xpers2 year,bw(8) lcol(blue*1.2)lpat(solid)col(blue*.25)legend(off) ///
			xtitle("Year", size(small)) saving(h1.gph,replace) ///
			ytitle(Marginal effect of security personalism)ylab(0 .02 .04 .06) ///
			yline(0,lcol(red))  xlab(1950(10)2010)tit(Onset years only)
			
		* All campaign years *
		tab any_repression hipers if xongoing==1,col
		krls any_repression xpers2 if xongoing==1  
		krls any_repression xpers2 lreduration lnregion lpop lnmembers lt period* if xongoing==1,d(k2)  
		krls any_repression xpers2 lreduration lnregion lpop lnmembers lt period* lag1repress lag2repress lag3repress if xongoing==1
		twoway lpolyci k2_xpers2 lreduration,bw(.5) lcol(blue*1.2)lpat(solid)col(blue*.25)legend(off) ///
			xtitle("Protest campaign duration (years, log scale)", size(small)) saving(h2.gph,replace) ///
			ytitle(Marginal effect of security personalism)ylab(0 .02 .04 .06) tit(All campaign years) ///
			yline(0,lcol(red))  xlab(0 "1" .69 "2"  1.39 "4"    2.079 "8" 2.48 "12")
		gr combine h1.gph h2.gph,xsize(8)
		graph export "$dir\navco-repression.pdf", as(pdf)   replace
		
	******************************************************
	***** Security forces defection during campaign ******
	******************************************************
		use temp-fe,clear
		xtset gwf_caseid year
		qui sum lnmembers if any_sec_def~=.
		gen xmembers = (lnmembers-r(mean))/ r(sd)
		qui sum xpers2
		gen ipers2 = (xpers2+abs(r(min)))/(r(max)+abs(r(min)))
		global x1 = "lt period2-period12 lpop xmembers lnregion repress lreduration"
		xi:logit any_sec_defect  xpers2 period2-period12,  vce(cluster id)
		est store def1
		xi:logit any_sec_defect  xpers2 lt period2-period12,  vce(cluster id)
		est store def2
		xi:logit any_sec_defect  xpers2 $x1,  vce(cluster id)
		est store def3
		xi:logit any_sec_defect  xpers2 $x1 milethnic_homo,  vce(cluster id)
		margins,dydx(xpers2 xmembers milethnic_homo) /* lnmembers & xpers2 both now Stdev==1 scale for margins */
		est store def4
		xi:logit any_sec_defect  ipers2 $x1 milethnic_homo,  vce(cluster id)
		sum ipers2 xmembers milethnic_homo if e(sample)==1
		margins,dydx(ipers2 xmembers milethnic_homo)  /* milethnic_homo & xpers2 both now (0,1) scale for margins */
		krls any_sec_defect xpers2 
		xi:krls any_sec_defect xpers2 $x1
		local c="lag_xongoing election civwar priordem coup coupA coupS excluded milethnic_het nmc_logmilper nmc_logmilex logoil loggdp xpers1 support"
		foreach v of local c {
			di "`v'"
			qui xi:logit any_sec_defect xpers2 $x1 i.region `v', cluster(id)
			lincom xpers2
			qui xi:krls any_sec_defect xpers2 $x1 i.region `v'
			matrix b = e(b)
			di b[1,1]
		}
		  * modeling unit hetero increases variance but also increases estimate size; so NOT modeling unit effect yields LOWER estimates, which are reported *
		xi:logit any_sec_defect xpers2 $x1,cluster(gwf_caseid)
		xi:xtlogit any_sec_defect xpers2 $x1,vce(cluster gwf_caseid)
		xi:xtlogit any_sec_defect xpers2 lt coldwar lpop xmembers lnregion repress lreduration,fe
		estout def1 def2 def3 def4 using TableE1.tex, cells(b(star  fmt(%9.3f)) se(par fmt(%9.3f))) ///
		stats(r2 N N_clust) style(tex) replace label starlevels(* 0.05) title(\label{tabE1})
		
		**** Alternative measures of personalization ****
		xi:logit any_sec_defect gwf_pers $x1,vce(cluster id)
		est store def_p1
		xi:logit any_sec_defect xpers $x1,vce(cluster id)
		est store def_p2
		xi:logit any_sec_defect xpers1 $x1,vce(cluster id)
		est store def_p3
		xi:logit any_sec_defect xpers2 $x1,vce(cluster id)
		est store def_p4
		xi:logit any_sec_defect xpers1 xpers2 $x1,vce(cluster id)
		est store def_p5


	********************************************************************** 
	******************     Democratization outcomes     ******************
	**********************************************************************
	use temp-fe,clear
	tab gdem gwf_case_fail
	forval i =1/12{
		egen m_period`i'=mean(period`i'),by(gwf_caseid)
	}
 	
	* show within transform equivalency *
	qui reghdfe gdem ld lnregion xpers2,a(gwf_caseid  year)cluster(gwf_caseid)
	lincom xpers2
	est store dem1
	xtsum gdem if e(sample)==1
	qui reghdfe gdem ld lnregion period* xpers2,a(gwf_caseid)cluster(gwf_caseid)
	lincom xpers2
	qui xtreg gdem ld lnregion period* xpers2,fe i(gwf_caseid)cluster(gwf_caseid)
	lincom xpers2
	qui reg gdem ld lnregion period* xpers2 m_ld m_lnregion m_coldwar m_xpers2 m_gdem m_period*,cluster(gwf_caseid)
	lincom xpers2
	qui reg gdem ld lnregion period* xpers2 m_ld m_lnregion m_coldwar m_xpers2 m_period*,cluster(gwf_caseid)
	lincom xpers2	
 
	  * CRE * Mundlak and Chamberlain within transformation, See Wooldridge 2002
	qui meprobit gdem ld lnregion period* xpers2 m_ld m_lnregion m_xpers2 m_period* || gwf_caseid: ,vce(cluster gwf_caseid)
	margins,dydx(xpers2) 
	est store dem2
	xtsum gdem if e(sample)==1
	qui probit gdem ld lnregion period* xpers2 m_ld m_lnregion m_xpers2  m_period*,cluster(gwf_caseid)
	margins,dydx(xpers2)
	
	* Add lagged trend in democracy *
	reghdfe gdem l1v2x_pol l2v2x_pol ld lnregion xpers2,a(gwf_caseid year)cluster(gwf_caseid)
	est store dem3
	xtsum gdem if e(sample)==1
		* Check with 3 & 4 year lags *
	reghdfe gdem l1v2x_pol l2v2x_pol l3v2x_pol l4v2x_pol ld lnregion xpers2,a(gwf_caseid year)cluster(gwf_caseid)
	egen m_l1v2x_pol=mean(l1v2x_pol) if sample==1,by(gwf_caseid)
	egen m_l2v2x_pol=mean(l2v2x_pol) if sample==1,by(gwf_caseid)
	meprobit gdem l1v2x_pol l2v2x_pol ld lnregion period* xpers2 ///
		m_ld m_lnregion m_xpers2 m_l1v2x_pol m_l2v2x_pol m_period* ///
		|| gwf_caseid: ,vce(cluster gwf_caseid)
	xtsum gdem if e(sample)==1
	est store dem4
	margins,dydx(xpers2) 
	
	* Interactive fixed effects *
	regife gdem ld lnregion xpers2 ,a(gwf_caseid year) factor(gwf_caseid year,1) vce(cluster gwf_caseid)
	est store dem5
	
	* Nonlinear case-specific time trends *
	xi:reg gdem i.gwf_caseid*time i.gwf_caseid*time2 ld lnregion xpers2 ,  vce(cluster gwf_leaderid)
	est store dem6
	
	estout dem1 dem2 dem3 dem4 dem5 using Table2.tex, cells(b(star  fmt(%9.4f)) se(par fmt(%9.3f))) ///
		stats(r2 N N_clust) style(tex) replace label starlevels(* 0.05) title(\label{tab2})
 
	* Interactive fixed effects, within: m_ are panel unit means; y_ are year means; and mXy_ are their interaction *
	meprobit gdem ld lnregion  xpers2 ///
		m_ld m_lnregion m_xpers2 y_ld y_lnregion y_xpers2  ///
		mXy_ld mXy_lnregion mXy_xpers2 ||  ///
		gwf_caseid: ,vce(cluster gwf_caseid)

	 * Check with IRT-2PL instead of generalized SFP *
	 reghdfe gdem ld lnregion irtpers,a(gwf_caseid year)cluster(gwf_caseid)
	 reghdfe gdem l1v2x_pol l2v2x_pol ld lnregion irtpers,a(gwf_caseid year)cluster(gwf_caseid)

	* Kernel estimator *
	krls gdem ld lnregion coldwar xpers2 m_ld m_lnregion m_coldwar m_xpers2, ///
		d(k)lambda(2.4)ltolerance(4.5)sigma(8)
	replace k_xpers2=k_xpers2*100
	*use temp-krls,clear
	twoway lpolyci k_xpers2 year if year>=1949, yline(0,lcol(red)lpat(dash))legend(off)ylab(-2(1)0) ///
		xlab(1950(10)2010)xtit(Year) ytit(Marginal effect of Security personalism) ///
		tit(Security personalism and democratic transition) ///
		note("Kernel regression, within transformation to proxy for regime-case fixed effects",pos(6)size(vsmall))
	graph export "$dir\krls-dem.pdf", as(pdf) replace  
	save temp-krls,replace
	
	 
		
		
	*** Other ways of modeling the time trend ***
	use temp-fe,clear
	forval i =1/12{
		egen m_period`i'=mean(period`i'),by(gwf_caseid)
	}
	qui reghdfe gdem ld lnregion coldwar xpers2,a(gwf_caseid)cluster(gwf_caseid)
	est store timedem1
	lincom xpers2
	qui reghdfe gdem ld lnregion d19* d20 xpers2,a(gwf_caseid)cluster(gwf_caseid)
	est store timedem2
	lincom xpers2
	qui reghdfe gdem ld lnregion period* xpers2,a(gwf_caseid)cluster(gwf_caseid)
	est store timedem3
	lincom xpers2
	qui reghdfe gdem ld lnregion time xpers2,a(gwf_caseid)cluster(gwf_caseid)
	est store timedem4
	lincom xpers2
	qui reghdfe gdem ld lnregion time time2 xpers2,a(gwf_caseid)cluster(gwf_caseid)
	est store timedem5
	lincom xpers2
	xtset gwf_caseid year
	regife gdem ld lnregion xpers2,a(gwf_caseid year)factor(gwf_caseid year,1)vce(cluster gwf_caseid)
	est store timedem6
		  * Nonlinear models *
	qui xtprobit gdem ld lnregion xpers2 m_ld m_lnregion m_xpers2 y_ld y_lnregion y_xpers2,vce(cluster gwf_caseid)
	est store timedem7
	 local var = "d1960 d1970 d1980 d1990 d2000 time time2"
	 foreach v of local var {
		egen m_`v' =mean(`v') if sample==1,by(gwf_caseid)
	 }
	qui xtprobit gdem ld lnregion xpers2 m_ld m_lnregion m_xpers2 d19* d20 m_d19* m_d20,vce(cluster gwf_caseid)
	est store timedem8
	qui xtprobit gdem ld lnregion xpers2 m_ld m_lnregion m_xpers2 period*,vce(cluster gwf_caseid)
	est store timedem9
	qui xtprobit gdem ld lnregion xpers2 m_ld m_lnregion m_xpers2 time m_time,vce(cluster gwf_caseid)
	est store timedem10
	qui xtprobit gdem ld lnregion xpers2 m_ld m_lnregion m_xpers2 time time2 m_time m_time2,vce(cluster gwf_caseid)
	est store timedem11
	qui xtprobit gdem ld lnregion xpers2 m_ld m_lnregion m_xpers2 ///
		y_ld y_lnregion y_xpers2 mXy_ld mXy_lnregion mXy_xpers2,vce(cluster gwf_caseid) 
	est store timedem12
	
	label var xpers2 " "
	coefplot (timedem1, msym(d))(timedem2, msym(t))(timedem3, msym(oh))(timedem4, msym(plus))(timedem5, msym(P)) ///
			(timedem6, msym(d)),title("Fixed effects linear probability model", size(med))  relocate(xpers2 = 1)  ///
			drop(lnregion coldwar _cons ld m_* y_* mXy_* d19* d20* period* time*) vertical xtit(Model specification,size(small)height(5)) ///
			order(xpers2) yline(0) grid(glcolor(gs13)) mfcolor(white) ylabel(-.04(.02)0,labsize(vsmall)) ///
			 levels(95 90)   xlab(.64 `""Year" "effects""' .78 `""Decade" "effects""' .925 `""Period" "effects""' 1.07 `""Linear" "time trend""' ///
			1.21 `""Non-linear" "time trend""' 1.35 `""Interactive" "fixed effects""',labsize(small) )  ///
			ysize(3) xsize(3.5) xscale(range(0.62 1.35))saving(g1, replace) ytit("Marginal effect estimate", size(small) ///
			height(4))legend(off)note("90 (thick) and 95 (thin) percent confidence intervals",size(small)ring(1) pos(6))
		
			 
		gen n=_n
		gen beta=.
		gen hi=.
		gen lo=.
		gen hi90=.
		gen lo90=.
		local i=1
		forval v =7/12 {
			di "`v'"
			est restore timedem`v'	 
			margins,dydx(xpers2)predict(pu0)  
			matrix beta =r(b)  
			local b = beta[1,1]
			qui replace beta=`b' if n==`i'
			margins,dydx(xpers2)predict(pu0)post
			matrix var = r(V) 
			local se =var[1,1]
			qui replace hi = `b' + sqrt(`se')*1.96 if n==`i'
			qui replace lo = `b' - sqrt(`se')*1.96 if n==`i'
			qui replace hi90 = `b' + sqrt(`se')*1.65 if n==`i'
			qui replace lo90 = `b' - sqrt(`se')*1.65 if n==`i'
			local i = `i' +1
		}
		graph twoway (rspike hi lo n if n<=6,lcolor(blue*.5)lwidth(thin) ///
			xtitle("Model specification",size(small)height(6)) ///
			yline(0,lp(dash)) ytitle("Marginal effect estimate", ///
			size(small) height(5))ylab(-.04(.02)0,)title("Correlated random effects probit",size(med))) ///
			(rspike hi90 lo90 n if n<=6,lcolor(blue*1) lwidth(medthick)) ///
			(scatter beta n if  n<=6,msym(dot)lpattern(solid)mcolor(blue*1.5) ///
			xscale(range(0.75 6.25))legend(lab(1 "95 ci")lab(3 "Estimate") size(small)pos(6)col(3)ring(1)) ///
			legend(order(3 1))saving(g2.gph,replace) xlab(1 `""Year" "effects""' 2 `""Decade" "effects""' ///
			3 `""Period" "effects""' 4 `""Linear" "time trend""' ///
			5 `""Non-linear" "time trend""' 6 `""Interactive" "fixed effects""'))
		gr combine g1.gph g2.gph,xsize(6)ysize(3)
					graph export "$dir\time-dem.pdf", as(pdf) replace  

			
	* Alternative Fixed effects *
	qui reghdfe gdem ld lnregion xpers2,a(gwf_caseid year)cluster(gwf_caseid)
	est store fedem1
	lincom xpers2
	qui reghdfe gdem ld lnregion xpers2,a(cowcode year)cluster(gwf_caseid)
	est store fedem2
	lincom xpers2
	qui reghdfe gdem lt lnregion xpers2,a(gwf_leaderid year)cluster(gwf_leaderid)
	est store fedem3
	lincom xpers2
	qui xi:xtreg gdem i.year ld lnregion xpers2,i(gwf_caseid)cluster(gwf_caseid)
	est store fedem4
	lincom xpers2
	qui xi:xtreg gdem i.year ld lnregion xpers2,i(cowcode)cluster(cowcode)
	est store fedem5
	lincom xpers2
	qui xi:xtreg gdem i.year lt lnregion xpers2,i(gwf_leaderid)cluster(gwf_leaderid)
	est store fedem6
	lincom xpers2
	
	coefplot (fedem1, msym(d))(fedem2, msym(t))(fedem3, msym(oh))(fedem4, msym(plus))(fedem5, msym(P)) ///
		(fedem6, msym(d)),title("Linear probability models", size(med))   ///
		drop(lnregion ld lt _Iyear_* _cons)   xtit(Model specification,size(small)height(5)) ///
		order(xpers2) xline(0) grid(glcolor(gs13)) mfcolor(white) xlabel(-.04(.02)0,labsize(vsmall)) ///
		levels(95 90)  ysize(3) xsize(3.5)   xtit("Marginal effect estimate", size(small) ///
		height(4))legend(lab(3 "Regime FE")lab(6 "Country FE")lab(9 "Leader FE") ////
		lab(12 "Regime RE")lab(15 "Country RE")lab(18 "Leader RE") ) ///
		note("90 (thick) and 95 (thin) percent confidence intervals; year effects in all specifications",size(vsmall)ring(1) pos(6))
	graph export "$dir\dem-effects.pdf", as(pdf) replace  
	
 
 
 	* All regime collapse *
	gen fail = gwf_case_fail
	recode fail (1=0) if gdem==1
 	qui reghdfe fail ld lnregion xpers2,a(gwf_caseid  year)cluster(gwf_caseid)
	lincom xpers2
	est store fail1
	xtsum fail if e(sample)==1
	xtset gwf_caseid year
	qui xtprobit fail ld lnregion period* xpers2 m_ld m_lnregion m_xpers2 m_period*, vce(cluster gwf_caseid)
	margins,dydx(xpers2) 
	lincom xpers2
	est store fail2
	xtsum fail if e(sample)==1
	estout fail1 fail2 using TableD1.tex, cells(b(star  fmt(%9.4f)) se(par fmt(%9.3f))) ///
		stats(r2 N N_clust) style(tex) replace label starlevels(* 0.05) title(\label{tabD1})
  
 	*** Cox duration model ***
	use temp-id,clear
	stset gwf_case_duration, failure(gdem) id(gwf_caseid)
	stcoxkm if gwf_case_duration<=100, by(treat) obs1opts(lcol(red)lpat(solid)mcol(red)msym(P)c(J)) ///
		obs2opts(lcol(blue)lpat(solid)mcol(blue)msym(P)c(J))  ///
		pred1opts(lcol(red)lpat(dash)mcol(red)msym(Oh)c(J)) ///
		pred2opts(lcol(blue)lpat(dash)mcol(blue)msym(Oh)c(J)) ///
		legend(pos(1)ring(0))ties(efron)xtit("Regime duration, years") ///
		note(Binary treatment is Security personalism greater than the sample median value,size(vsmall)pos(6))
		graph export "$dir\dem-coxkm.pdf", as(pdf)   replace
	stcox period* lnregion xpers2,vce(cluster gwf_caseid)nohr
	est store demcox1
		* Cox, within transform *
	forval i =1/12{
		egen m_period`i'=mean(period`i'),by(gwf_caseid)
	}
	stcox lnregion period* xpers2 m_lnregion m_period* m_xpers2,vce(cluster gwf_caseid)	nohr
	est store demcox2
		* Cox, with shared frailty by regime *
	gen ldXlnregion=ld*lnregion
	stcox period* lnregion xpers2 ldXlnregion,   shared(gwf_caseid)nohr
	est store demcox3
	estat phtest, rank detail
		* Cox, with strata by country *
	stcox period* lnregion xpers2, vce(cluster gwf_caseid)  strata(cowcode)	nohr
	estat phtest, rank detail
	est store demcox4
	stcox period* lnregion lt treat, vce(cluster gwf_caseid)  strata(cowcode)nohr
	estat phtest, rank detail
	stcox period* lnregion xpers2 loggdp lpop logoil, vce(cluster gwf_caseid)  strata(cowcode)nohr
	est store demcox5
	estat phtest, rank detail
	stcox period* lnregion xpers2 loggdp lpop logoil l1v2x_polyarchy l2v2x_polyarchy, vce(cluster gwf_caseid)   strata(cowcode)nohr	
	est store demcox6
	estout demcox* using TableD2.tex, cells(b(star  fmt(%9.4f)) se(par fmt(%9.3f))) ///
		stats(r2 N N_clust) style(tex) replace label starlevels(* 0.05) title(\label{tabD2})
  	
 
		**** Adding covariates, one at a time ****
		use temp-fe,clear
		recode debruin_affcount (5 4 3=3)
		spearman xpers2 v2x_clpol v2clkill v2cltort v2juhcind v2x_jucon v2x_frassoc_thick v2x_freexp_altinf v2x_partipdem v2x_libdem v2x_polyarchy
		gen n=_n
		gen beta=.
		gen hi=.
		gen lo=.
		gen hi90=.
		gen lo90=.
		gen varname=""
		xtset gwf_caseid year
		qui reghdfe gdem ld lnregion xpers2,a(gwf_caseid  year)cluster(gwf_caseid)
		nlcom _b[xpers2],post
		matrix beta =e(b)  
		global b = beta[1,1]	
		global num =40
		global cvar ="ld lnregion"
		global varA= "lpop loggdp logoil support gwf_mil leadermil seizure_coup coup12 election ythbul4 wdipopurbmi wditrade"
		global varB= "polparcomp polpolcomp polexconst l1v2x_polyarchy l2v2x_polyarchy l1v2x_partipdem l1v2x_freexp_altinf l1v2x_frassoc_thick "
		global varC= "l1v2x_clpol l1v2x_jucon l1v2juhcind l1e_v2x_neopat priordem legcomp inst excluded monoethnic multiethnic civwar grow"
		global varD= "milethnic_homo nmc_logmilper nmc_logmilex effectivenumber debruin_cbcount debruin_ha_cbcount debruin_affcount lag_xongoing"
		local var = "$varA $varB $varC $varD"
		local i =1
		foreach v of local var {
			di "`v'"
			qui xtset gwf_caseid
			qui reghdfe gdem `v' ld lnregion xpers2,a(gwf_caseid  year)cluster(gwf_caseid)
			qui nlcom _b[xpers2],post
			matrix beta =e(b)  
			local b = beta[1,1]
			qui replace beta=`b' if n==`i'
			matrix var = e(V) 
			local se =var[1,1]
			qui replace hi = `b' + sqrt(`se')*1.96 if n==`i'
			qui replace lo = `b' - sqrt(`se')*1.96 if n==`i'
			qui replace hi90 = `b' + sqrt(`se')*1.65 if n==`i'
			qui replace lo90 = `b' - sqrt(`se')*1.65 if n==`i'
			qui replace varname = "`v'" if n==`i'
			local i = `i' +1
		}
		label define varlab 1 "Population" 2 "GDP per capita (log)" 3 "Oil per capita (log)"  ///
			4 "Support party" 5 "Military junta" 6 "Military leader" 7 "Coup seizure" 8 "Recent coup" 9 "Election" ///
			10 "Youth bulge" 11 "Urban pop." 12 "Trade"   13 "Parcomp" 14 "Polcomp" 15 "Exec. constr." ///
			16 "V-Polyarchy,t-1" 17 "V-Polyarchy,t-2" 18 "V-Participation,t-1" 19 "V-Free express. info,t-1" ///
			20 "V-Free express. org,t-1" 21 "V-Civil lib.,t-1" 22 "V-Judical constr.,t-1" 23 "V-Judicial indep.,t-1" ///
			24 "V-Neopatrimonialism,t-1" 25 "Prior democracy" ///
			26 "Leg. comp." 27 "Institutions" 28 "Excluded pop" 29 "Monoethnic party" 30 "Multiethnic party" ///
			31 "Civil war" 32 "Growth" 33 "Ethnic homo military"  34 "Military personnel (log)" ///
			35 "Military exp. (log)" 36 "Effective number" 37 "Counter-weights" 38 "Heavily armed counter-weights" ///
			39 "Affiliated paramilitaries" 40 "Ongoing protests,t-1",replace
		label values n varlab
		twoway (scatter beta n if n<=$num,mcol(blue)yline($b,lcol(red)lpat(dash_dot))) ///
			(rspike hi lo n if n<=$num,lw(vthin)lcol(blue)) ///
			(rspike hi90 lo90 n if n<=$num,lcol(blue)lw(medium)ytitle("{&beta}{sub:Security personalization}", ///
			size(large)height(4)) ylab(-.06(.02)0) ///
			xtitle(Added covariate)yline(0,lpat(dash)lcol(gs6))xlab(1(1)$num,valuelabel angle(90))legend(off))
		graph export "$dir\dem-baseline-covariates.pdf", as(pdf)   replace
		
		*** IV-2SLS ***
		use temp-fe,clear
		xtset gwf_caseid year
		forval i =1/4 {
			gen l`i'vdem = l`i'.v2x_polyarchy
		}
		gen militrank2=militrank^2
		xi:ivreghdfe gdem ld lnregion (xpers2=militrank militrank2),absorb(gwf_caseid year)cluster(gwf_caseid) 
		est store demiv1
		xi:ivreghdfe gdem ld lnregion l1vdem l2vdem (xpers2=militrank militrank2),absorb(gwf_caseid year)cluster(gwf_caseid)   
		est store demiv2
		
		xi:ivreg2h gdem i.year  ld lnregion (xpers2=militrank militrank2),fe cluster(gwf_caseid) partial(i.year) gmm2s
		est store demiv3
		xi:ivreg2h gdem i.year l1vdem l2vdem  ld lnregion (xpers2=militrank militrank2),fe cluster(gwf_caseid) partial(i.year) gmm2s
		est store demiv4
		
		* Check with 3- and 4-lags *
		xi:ivreg2h gdem i.year l1vdem l2vdem l3vdem  ld lnregion (xpers2=militrank militrank2),fe cluster(gwf_caseid) partial(i.year) gmm2s
		xi:ivreg2h gdem i.year l1vdem l2vdem l3vdem l4vdem  ld lnregion (xpers2=militrank militrank2),fe cluster(gwf_caseid) partial(i.year) gmm2s
		
		estout demiv* using TableD3.tex, cells(b(star  fmt(%9.4f)) se(par fmt(%9.3f))) ///
			stats(r2 N N_clust widstat) style(tex) replace label starlevels(* 0.05) title(\label{tabD3})
			
			
		xi:reghdfe gdem lnregion l1vdem l2vdem l3vdem l4vdem xpers2,absorb(gwf_case_duration gwf_caseid year)cluster(gwf_caseid)
		xi:reghdfe gdem lnregion lpop loggdp logoil l1vdem l2vdem l3vdem l4vdem l.xpers2 xpers2,absorb(gwf_case_duration gwf_caseid year)cluster(gwf_caseid)
		xi:reghdfe gdem lnregion lpop loggdp logoil l1vdem l2vdem l3vdem l4vdem xpers2,absorb(gwf_case_duration gwf_caseid year)cluster(gwf_caseid)
	
	
	********* VDEM democracy score ************
	  reghdfe v2x_polyarchy ld lnregion xpers2,a(gwf_caseid year)cluster(gwf_leaderid)
	  egen minyr= min(year),by(gwf_caseid)
	  gen ovdem = l1v2x_poly if year==minyr
	  egen ivdem = max(ovdem),by(gwf_caseid)
	  reghdfe v2x_polyarchy ivdem ld lnregion xpers2,a(cowcode year)cluster(gwf_leaderid)

log close
		
 
		  
		****************** THE END *****************
		
